set.seed(1)
# Data manipulation & cleaning
library(tidyverse)
library(reshape2)
library(nlme)        # gapply
library(lubridate)   # date cleaning

# Viz
library(gridExtra)
library(ggmap)
library(maps)
library(mapdata)

# Modeling & stats
library(mgcv)       # gam
library(R0)
library(Kendall)    # Mann-Kendall

# Time series & forecasting
library(forecast)
library(zoo)
library(astsa)
library(fpp2)

TODO

RT ANALYSIS

#read in full DSHS data with county and TSA information
# TOOD: rename
full.dshs <- read.csv("combined-datasets/county.csv")
tsa.dshs <- read.csv("combined-datasets/tsa.csv")
phr.df <- read.csv("combined-datasets/phr.csv")
rt.data.clean<-function(covid.data) {
  
  # select relevant metadata & DSHS population estimates
  rt.cols <- c("County", "Date", "TSA_Name", "TSA", "PHR", "PHR_Name",
               "Metro_Area", "Cases_Daily", "Population_DSHS")
  
  covid.data <- covid.data[, names(covid.data) %in% rt.cols]
             
  # set correct types
  covid.data$Date <- as.Date(covid.data$Date)
  covid.data$Cases_Daily <- as.numeric(covid.data$Cases_Daily)
  
  # TODO: decide on how to handle drops in cumulative counts -> negative daily cases
  covid.data <- covid.data %>% 
    mutate(Cases_Daily = ifelse(Cases_Daily < 0, 0, Cases_Daily))
  
  return(covid.data)
}
#separate county, metro, TSA and state data into separate dataframes
rt.data.organize <- function(mycounty, mytsa, myphr){
  
  ### COUNTY ###
  cleaned.county <- rt.data.clean(mycounty)
  county.df <- cleaned.county%>%dplyr::select(-TSA_Name, -Metro_Area, -TSA)
  
  ### TSA ###
  TSA.df <- rt.data.clean(mytsa)

  ### PHR ###
  PHR.df <- rt.data.clean(myphr)
  ### METRO ###
  # calculate daily cases and population at metro level, drop repeated rows
  metro.temp <- cleaned.county %>%
    group_by(Metro_Area,Date) %>%
    mutate(metro_cases_daily=sum(Cases_Daily, na.rm=TRUE)) %>%
    mutate(metro_pop_DSHS=sum(Population_DSHS, na.rm=TRUE)) %>%
    dplyr::select(Date, Metro_Area, metro_cases_daily, metro_pop_DSHS) %>%
    distinct()
  
  
  metro.df <- data.frame(Date = metro.temp$Date,
                         Metro_Area = metro.temp$Metro_Area,
                         Cases_Daily = metro.temp$metro_cases_daily,
                         Population_DSHS = metro.temp$metro_pop_DSHS)
  
  
  ### STATE ###
  # calculate daily cases and population at state level, drop repeated rows
  state.temp <- TSA.df %>%
    group_by(Date) %>%
    mutate(state_cases_daily=sum(Cases_Daily, na.rm=TRUE)) %>%
    mutate(state_pop_DSHS=sum(Population_DSHS, na.rm=TRUE)) %>%
    dplyr::select(Date, state_cases_daily, state_pop_DSHS) %>%
    distinct()
  
  
  state.df <- data.frame(Date = state.temp$Date, 
                         Cases_Daily = state.temp$state_cases_daily,
                         Population_DSHS = state.temp$state_pop_DSHS)

  ### OUTPUT ###
  rt.df.out <- list(county=county.df, TSA=TSA.df, state=state.df, metro=metro.df, phr = PHR.df)
  return(rt.df.out)
}
rt.df.extraction <- function(Rt.estimate.output) {

  # extract r0 estimate values into dataframe
  rt.df <- setNames(stack(Rt.estimate.output$estimates$TD$R)[2:1], c('Date', 'Rt'))
  rt.df$Date <- as.Date(rt.df$Date)

  # get 95% CI
  CI.lower.list <- Rt.estimate.output$estimates$TD$conf.int$lower
  CI.upper.list <- Rt.estimate.output$estimates$TD$conf.int$upper

  #use unlist function to format as vector
  CI.lower <- unlist(CI.lower.list, recursive = TRUE, use.names = TRUE)
  CI.upper <- unlist(CI.upper.list, recursive = TRUE, use.names = TRUE)

  rt.df$lower <- CI.lower
  rt.df$upper <- CI.upper
  
  rt.df <- rt.df %>%
    mutate(lower = replace(lower, Rt == 0, NA)) %>%
    mutate(upper = replace(upper, Rt == 0, NA)) %>%
    mutate(Rt = replace(Rt, Rt == 0, NA))
  
  return(rt.df)
}
#Function to generate Rt estimates, need to read in data frame and population size
covid.rt <- function(mydata) {
  
  ### DECLARE VALS ###
  
  #set generation time
  #Tapiwa, Ganyani "Esimating the gen interval for Covid-19":

  # LOGNORMAL OPS
  # gen.time<-generation.time("lognormal", c(4.0, 2.9))
  # gen.time<-generation.time("lognormal", c(4.7,2.9)) #Nishiura

  # GAMMA OPS
  # gen.time<-generation.time("gamma", c(5.2, 1.72)) #Singapore
  # gen.time<-generation.time("gamma", c(3.95, 1.51)) #Tianjin
  gen.time <- generation.time("gamma", c(3.96, 4.75))
  print(as.character(mydata[1,2]))
  
  # TODO: consider removing this and handling na cases directly
  #change na values to 0
  mydata <- mydata %>% replace(is.na(.),0)
  sum.daily.cases <- sum(mydata$Cases_Daily)
  pop.DSHS <- mydata$Population_DSHS[1]
  
  #Get 7 day moving average of daily cases
  mydata$MA_7day<-rollmean(mydata$Cases_Daily, k=7, na.pad=TRUE, align='right')
  #change na values to 0
  mydata<-mydata%>%replace(is.na(.),0)

  #create a vector of new cases 7 day moving average
  mydata.new<-pull(mydata, MA_7day)
  #mydata.new <- pull(mydata, Cases_Daily)
  
  # get dates as vectors
  date.vector <- pull(mydata, Date)
  
  #create a named numerical vector using the date.vector as names of new cases
  #Note: this is needed to run R0 package function estimate.R()
  names(mydata.new) <- c(date.vector)
  
  #Find min.date, max.date and row number of max.date
  min.date <- min(mydata$Date)
  max.date <- max(mydata$Date)

  #get row number of March 15 and first nonzero entry
  #NOTE: for 7 day moving average start March 15, for daily start March 9
  #find max row between the two (this will be beginning of rt data used)
  march15.row <- which(mydata$Date=="2020-03-15")
  first.nonzero <- min(which(mydata$Cases_Daily>0))
  last.nonzero <- max(which(mydata$Cases_Daily>0))
  minrow <- max(march15.row, first.nonzero)

  ### R0 ESTIMATION ###
  # TODO: establish better criteria for # of required daily cases
  # mean daily cases; average last x number of days; % of region pop etc.
  if(!is.na(minrow) & sum.daily.cases > 100) {
    
    # declare limits for rt estimation (first nonzero date & last nonzero date)
    i <- minrow
    j <- as.integer(min(last.nonzero, max.date))

    #reduce the vector to be rows from min date (March 9 or first nonzero case) to current date
    mydata.newest <- mydata.new[i:j]
    rt.DSHS <- estimate.R(mydata.newest, 
                        gen.time, 
                        begin=as.integer(1),
                        end=length(mydata.newest),
                        methods=c("TD"), 
                        pop.size=pop.DSHS,
                        nsim=1000)
  
    rt.DSHS.df <- rt.df.extraction(rt.DSHS)
    
  } else {
    # error catch small regions & na values for minrow
    rt.DSHS.df <- data.frame(Date = as.Date(mydata$Date),
                             Rt = rep(NA, length(mydata$Date)),
                             lower = rep(NA, length(mydata$Date)),
                             upper = rep(NA, length(mydata$Date)))
  }
  
  return(rt.DSHS.df)
}

County

#apply the covid.rt function to the county.datily data by county 
county.rt.output <- nlme::gapply(county.daily, FUN = covid.rt, groups = county.daily$County)
## [1] "Anderson"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Andrews"
## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.142857142857143, `2020-04-05` =
## 0.857142857142857, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-04` = 0.142857142857143, `2020-04-05` =
## 0.857142857142857, : Using initial incidence as initial number of cases.
## [1] "Angelina"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Aransas"
## [1] "Archer"
## [1] "Armstrong"
## [1] "Atascosa"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Austin"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Bailey"
## Warning in est.R0.TD(epid = c(`2020-05-05` = 0.142857142857143, `2020-05-06` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-05-05` = 0.142857142857143, `2020-05-06` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Bandera"
## [1] "Bastrop"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Baylor"
## [1] "Bee"
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Bell"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Bexar"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Blanco"
## [1] "Borden"
## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf
## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf
## [1] "Bosque"
## [1] "Bowie"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Brazoria"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Brazos"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Brewster"
## Warning in est.R0.TD(epid = c(`2020-05-01` = 0.142857142857143, `2020-05-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-05-01` = 0.142857142857143, `2020-05-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Briscoe"
## [1] "Brooks"
## [1] "Brown"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Burleson"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Burnet"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Caldwell"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Calhoun"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Callahan"
## [1] "Cameron"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Camp"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Carson"
## [1] "Cass"
## [1] "Castro"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Chambers"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Cherokee"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Childress"
## [1] "Clay"
## [1] "Cochran"
## [1] "Coke"
## [1] "Coleman"
## [1] "Collin"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Using initial incidence as initial number of cases.
## [1] "Collingsworth"
## [1] "Colorado"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.285714285714286, `2020-04-02` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.285714285714286, `2020-04-02` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Comal"
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.428571428571429, `2020-03-23` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.428571428571429, `2020-03-23` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Comanche"
## [1] "Concho"
## [1] "Cooke"
## Warning in est.R0.TD(epid = c(`2020-04-10` = 0.142857142857143, `2020-04-11` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-10` = 0.142857142857143, `2020-04-11` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Coryell"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Cottle"
## [1] "Crane"
## [1] "Crockett"
## [1] "Crosby"
## [1] "Culberson"
## [1] "Dallam"
## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Dallas"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 1.14285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 1.14285714285714, : Using initial incidence as initial number of cases.
## [1] "Dawson"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Deaf Smith"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Delta"
## [1] "Denton"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "DeWitt"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.142857142857143, `2020-03-20` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.142857142857143, `2020-03-20` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Dickens"
## [1] "Dimmit"
## [1] "Donley"
## [1] "Duval"
## [1] "Eastland"
## [1] "Ector"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Edwards"
## [1] "El Paso"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Ellis"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Erath"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Falls"
## [1] "Fannin"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Fayette"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Fisher"
## [1] "Floyd"
## [1] "Foard"
## [1] "Fort Bend"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 1.14285714285714, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Franklin"
## [1] "Freestone"
## [1] "Frio"
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Gaines"
## [1] "Galveston"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Garza"
## [1] "Gillespie"
## [1] "Glasscock"
## [1] "Goliad"
## [1] "Gonzales"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Gray"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Grayson"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Gregg"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Grimes"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Guadalupe"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.571428571428571, `2020-03-26` =
## 1.14285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.571428571428571, `2020-03-26` =
## 1.14285714285714, : Using initial incidence as initial number of cases.
## [1] "Hale"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hall"
## [1] "Hamilton"
## [1] "Hansford"
## [1] "Hardeman"
## [1] "Hardin"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Harris"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.857142857142857, `2020-03-16` =
## 0.714285714285714, : Using initial incidence as initial number of cases.
## [1] "Harrison"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hartley"
## [1] "Haskell"
## [1] "Hays"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hemphill"
## [1] "Henderson"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hidalgo"
## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.285714285714286, `2020-03-24` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.285714285714286, `2020-03-24` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Hill"
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hockley"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Hood"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Hopkins"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Houston"
## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Howard"
## [1] "Hudspeth"
## [1] "Hunt"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Hutchinson"
## [1] "Irion"
## [1] "Jack"
## [1] "Jackson"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Jasper"
## [1] "Jeff Davis"
## [1] "Jefferson"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 1, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 1, : Using initial incidence as initial number of cases.
## [1] "Jim Hogg"
## [1] "Jim Wells"
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Johnson"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Jones"
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-08` = 0.142857142857143, `2020-04-09` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Karnes"
## [1] "Kaufman"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.285714285714286, `2020-03-25` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Kendall"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Kenedy"
## [1] "Kent"
## [1] "Kerr"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Kimble"
## [1] "King"
## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf
## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf
## [1] "Kinney"
## [1] "Kleberg"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Knox"
## [1] "La Salle"
## [1] "Lamar"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Lamb"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.142857142857143, `2020-03-29` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Lampasas"
## [1] "Lavaca"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Lee"
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.142857142857143, `2020-04-03` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Leon"
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Liberty"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Limestone"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Lipscomb"
## [1] "Live Oak"
## [1] "Llano"
## [1] "Loving"
## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf
## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf
## [1] "Lubbock"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Lynn"
## [1] "Madison"
## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-15` = 0.142857142857143, `2020-04-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Marion"
## [1] "Martin"
## [1] "Mason"
## [1] "Matagorda"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Maverick"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "McCulloch"
## [1] "McLennan"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Using initial incidence as initial number of cases.
## [1] "McMullen"
## [1] "Medina"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Menard"
## [1] "Midland"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Milam"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.285714285714286, `2020-03-26` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Mills"
## [1] "Mitchell"
## [1] "Montague"
## [1] "Montgomery"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Moore"
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Morris"
## [1] "Motley"
## [1] "Nacogdoches"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Navarro"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Newton"
## [1] "Nolan"
## [1] "Nueces"
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Ochiltree"
## [1] "Oldham"
## [1] "Orange"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Palo Pinto"
## [1] "Panola"
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.428571428571429, `2020-04-03` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-02` = 0.428571428571429, `2020-04-03` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "Parker"
## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.142857142857143, `2020-03-24` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-23` = 0.142857142857143, `2020-03-24` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Parmer"
## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-19` = 0.142857142857143, `2020-04-20` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Pecos"
## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-06` = 0.142857142857143, `2020-04-07` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Polk"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Potter"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.285714285714286, `2020-03-22` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.285714285714286, `2020-03-22` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Presidio"
## [1] "Rains"
## [1] "Randall"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.428571428571429, `2020-03-29` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.428571428571429, `2020-03-29` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Reagan"
## [1] "Real"
## [1] "Red River"
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Reeves"
## [1] "Refugio"
## [1] "Roberts"
## [1] "Robertson"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Rockwall"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Runnels"
## [1] "Rusk"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Sabine"
## [1] "San Augustine"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "San Jacinto"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.142857142857143, `2020-03-30` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "San Patricio"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "San Saba"
## [1] "Schleicher"
## [1] "Scurry"
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-16` = 0.142857142857143, `2020-04-17` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Shackelford"
## [1] "Shelby"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Sherman"
## [1] "Smith"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Somervell"
## [1] "Starr"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.285714285714286, `2020-03-28` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Stephens"
## [1] "Sterling"
## [1] "Stonewall"
## Warning in min(which(mydata$Cases_Daily > 0)): no non-missing arguments to min;
## returning Inf
## Warning in max(which(mydata$Cases_Daily > 0)): no non-missing arguments to max;
## returning -Inf
## [1] "Sutton"
## [1] "Swisher"
## [1] "Tarrant"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Taylor"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Terrell"
## [1] "Terry"
## [1] "Throckmorton"
## [1] "Titus"
## Warning in est.R0.TD(epid = c(`2020-04-03` = 0.142857142857143, `2020-04-04` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-03` = 0.142857142857143, `2020-04-04` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Tom Green"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Travis"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Trinity"
## [1] "Tyler"
## [1] "Upshur"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Upton"
## [1] "Uvalde"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Val Verde"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.142857142857143, `2020-03-27` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "Van Zandt"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Victoria"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.428571428571429, `2020-03-27` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.428571428571429, `2020-03-27` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Walker"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Waller"
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.428571428571429, `2020-03-30` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-29` = 0.428571428571429, `2020-03-30` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Ward"
## [1] "Washington"
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.714285714285714, `2020-03-29` =
## 0.714285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-28` = 0.714285714285714, `2020-03-29` =
## 0.714285714285714, : Using initial incidence as initial number of cases.
## [1] "Webb"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Wharton"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.428571428571429, `2020-03-26` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Wheeler"
## [1] "Wichita"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Wilbarger"
## [1] "Willacy"
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-27` = 0.142857142857143, `2020-03-28` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Williamson"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 1.14285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.285714285714286, `2020-03-20` =
## 1.14285714285714, : Using initial incidence as initial number of cases.
## [1] "Wilson"
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-25` = 0.142857142857143, `2020-03-26` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Winkler"
## [1] "Wise"
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-30` = 0.142857142857143, `2020-03-31` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Wood"
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-01` = 0.142857142857143, `2020-04-02` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Yoakum"
## [1] "Young"
## [1] "Zapata"
## Warning in est.R0.TD(epid = c(`2020-04-07` = 0.142857142857143, `2020-04-08` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-04-07` = 0.142857142857143, `2020-04-08` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "Zavala"
#convert list of dataframes to single dataframe with county as id col
county.rt.df <- data.table::rbindlist(county.rt.output, idcol = TRUE)
colnames(county.rt.df)[1] = 'County'

Metro

#apply the covid.rt function to the metro.daily data by metro region
metro.rt.output <- nlme::gapply(metro.daily, FUN=covid.rt, groups=metro.daily$Metro_Area)
## [1] "Metro"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 7, `2020-03-16` = 6, `2020-03-17` =
## 6.57142857142857, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 7, `2020-03-16` = 6, `2020-03-17` =
## 6.57142857142857, : Using initial incidence as initial number of cases.
## [1] "Non-Metro"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.285714285714286, `2020-03-16` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
#convert list of data frames to single data frame with metro as id col
metro.rt.df <- data.table::rbindlist(metro.rt.output, idcol=TRUE)
colnames(metro.rt.df)[1] = 'Metro_Area'

TSA

# plot with shaded confidence intervals
rt.plot<-function(rt.data, caption){
  library(ggplot2)
    caption<-rt.data$TSA[1]
  plot<-ggplot(rt.data, aes(Date, Rt))+
    geom_ribbon(aes(ymin=lower, ymax=upper), fill="light grey")+
    geom_point(color="blue", size=0.2)+
    labs(title=caption)
  plot
}
#apply the covid.rt function to the TSA.daily data by TSA region
TSA.rt.output <- nlme::gapply(TSA.daily, FUN=covid.rt, groups=TSA.daily$TSA)
## [1] "A"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.571428571428571, `2020-03-22` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.571428571428571, `2020-03-22` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "B"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.714285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 0.714285714285714, : Using initial incidence as initial number of cases.
## [1] "C"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "D"
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-21` = 0.142857142857143, `2020-03-22` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "E"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Using initial incidence as initial number of cases.
## [1] "F"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "G"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "H"
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-26` = 0.285714285714286, `2020-03-27` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "I"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "J"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.285714285714286, : Using initial incidence as initial number of cases.
## [1] "K"
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-24` = 0.142857142857143, `2020-03-25` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "L"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "M"
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-19` = 0.857142857142857, `2020-03-20` =
## 1, : Using initial incidence as initial number of cases.
## [1] "N"
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-18` = 0.142857142857143, `2020-03-19` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "O"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "P"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "Q"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.57142857142857, `2020-03-16` =
## 1.71428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.57142857142857, `2020-03-16` =
## 1.71428571428571, : Using initial incidence as initial number of cases.
## [1] "R"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.428571428571429, `2020-03-16` =
## 0.428571428571429, : Using initial incidence as initial number of cases.
## [1] "S"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "T"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "U"
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-22` = 0.142857142857143, `2020-03-23` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "V"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.142857142857143, `2020-03-21` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
#convert list of data frames to single data frame with TSA as variable
TSA.rt.df <- data.table::rbindlist(TSA.rt.output, idcol=TRUE)

# add TSA_Name & fix dates
colnames(TSA.rt.df)[1] <- 'TSA'
TSA.rt.df$Date <- as.Date(TSA.rt.df$Date)
TSA.rt.df <- merge(TSA.rt.df, TSA.daily[, c('TSA', 'TSA_Name', 'Date')], by = c('TSA', 'Date'))

# combine TSA 
TSA.rt.df$TSA <- paste0(TSA.rt.df$TSA, ' - ', TSA.rt.df$TSA_Name)
TSA.rt.df$TSA_Name <- NULL

#plot Rt by TSA
nlme::gapply(TSA.rt.df, FUN=rt.plot, groups=TSA.rt.df$TSA)
## $`A - Amarillo`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`B - Lubbock`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`C - Wichita Falls`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`D - Abilene`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`E - Dallas/Ft. Worth`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`F - Paris`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`G - Longview/Tyler`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`H - Lufkin`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`I - El Paso`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`J - Midland/Odessa`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`K - San Angelo`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`L - Belton/Killeen`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`M - Waco`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`N - Bryan/College Station`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`O - Austin`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`P - San Antonio`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`Q - Houston`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`R - Galveston`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`S - Victoria`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`T - Laredo`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`U - Corpus Christi`
## Warning: Removed 1 rows containing missing values (geom_point).

## 
## $`V - Lower Rio Grande Valley`
## Warning: Removed 1 rows containing missing values (geom_point).

PHR

phr.rt.output <- nlme::gapply(phr.daily, FUN=covid.rt, groups=phr.daily$PHR)
## [1] "1"
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 1.28571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-20` = 0.285714285714286, `2020-03-21` =
## 1.28571428571429, : Using initial incidence as initial number of cases.
## [1] "11"
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-17` = 0.142857142857143, `2020-03-18` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
## [1] "2/3"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 2.42857142857143, `2020-03-16` =
## 2.28571428571429, : Using initial incidence as initial number of cases.
## [1] "4/5N"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "6/5S"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 3, `2020-03-16` =
## 2.14285714285714, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 3, `2020-03-16` =
## 2.14285714285714, : Using initial incidence as initial number of cases.
## [1] "7"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "8"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.571428571428571, `2020-03-16` =
## 0.571428571428571, : Using initial incidence as initial number of cases.
## [1] "9/10"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 0.142857142857143, `2020-03-16` =
## 0.142857142857143, : Using initial incidence as initial number of cases.
#convert list of data frames to single data frame with phr as id col
phr.rt.df <- data.table::rbindlist(phr.rt.output, idcol=TRUE)

# add PHR_Name & fix dates
colnames(phr.rt.df)[1] <- 'PHR'
phr.rt.df$Date <- as.Date(phr.rt.df$Date)
phr.rt.df <- merge(phr.rt.df, phr.daily[, c('PHR', 'PHR_Name', 'Date')], by = c('PHR', 'Date'))

# combine phr 
phr.rt.df$PHR <- paste0(phr.rt.df$PHR, ' - ', phr.rt.df$PHR_Name)
phr.rt.df$PHR_Name <- NULL

State

state.rt.df <- covid.rt(state.daily)
## [1] "0"
## Warning in est.R0.TD(epid = c(`2020-03-15` = 7.28571428571429, `2020-03-16` =
## 6.28571428571429, : Simulations may take several minutes.
## Warning in est.R0.TD(epid = c(`2020-03-15` = 7.28571428571429, `2020-03-16` =
## 6.28571428571429, : Using initial incidence as initial number of cases.

export

write.csv(county.rt.df,"statistical-output/rt/county_rt.csv", row.names = FALSE)
write.csv(metro.rt.df, "statistical-output/rt/metro_rt.csv", row.names = FALSE)
write.csv(TSA.rt.df, "statistical-output/rt/tsa_rt.csv", row.names = FALSE)
write.csv(phr.rt.df, "statistical-output/rt/phr_rt.csv", row.names = FALSE)
write.csv(state.rt.df, "statistical-output/rt/state_rt.csv", row.names = FALSE)

grouping

colnames(county.rt.df)[1] <- 'Level'
colnames(metro.rt.df)[1] <- 'Level'
colnames(TSA.rt.df)[1] <- 'Level'
colnames(phr.rt.df)[1] <- 'Level'
state.rt.df$Level <- 'Texas'

county.rt.df$Level_Type = 'County'
metro.rt.df$Level_Type = 'Metro'
TSA.rt.df$Level_Type = 'TSA'
phr.rt.df$Level_Type = 'PHR'
state.rt.df$Level_Type = 'State'

combined.rt.df <- rbind(subset(county.rt.df, Date != max(Date)),
                        subset(metro.rt.df, Date != max(Date)),
                        subset(TSA.rt.df, Date != max(Date)),
                        subset(state.rt.df, Date != max(Date)),
                        subset(phr.rt.df, Date != max(Date)))

write.csv(combined.rt.df, 'statistical-output/rt/stacked_rt.csv', row.names = FALSE)
write.csv(combined.rt.df, 'tableau/stacked_rt.csv', row.names = FALSE)

TIME SERIES

#ts.rt.data.clean function will clean data, dropping unwanted variables
ts.data.clean<-function(covid.data) {

  ts_cols <- c("County", "Date", "TSA", "TSA_Name", "PHR", "PHR_Name",
               "Metro_Area", "Cases_Daily", "Cases_Cumulative",
               "Population_DSHS", "Hospitalizations_Total")
  
  covid.data <- covid.data[, names(covid.data) %in% ts_cols]
  
  # set correct types
  covid.data$Date <- as.Date(covid.data$Date)
  covid.data$Cases_Daily <- as.numeric(covid.data$Cases_Daily)
  
  #change any Cases_Daily below zero to zero
  covid.data <- covid.data %>% mutate(Cases_Daily = ifelse(Cases_Daily < 0, 0, Cases_Daily))
  return(covid.data)
}
#separate county, metro, TSA and State data into separate Data frames
ts.data.organize<-function(mycounty, mytsa, myphr){

  ### COUNTY ###
  # clean county vals and restrict to first date of collection
  cleaned.county <- ts.data.clean(mycounty)
  cleaned.county <- subset(cleaned.county, !is.na(Date) & Date >= as.Date('2020-03-04'))

  # Set column names
  county.df <- cleaned.county%>%dplyr::select(-TSA_Name, -Metro_Area)
  county.df$Level <- 'County'
  colnames(county.df)[2] = 'Level_Name'
  
  #change hospitalizations to numeric
  mytsa$Hospitalizations_Total<-as.numeric(mytsa$Hospitalizations_Total)
  
  
  # PHR 
  PHR.df <- ts.data.clean(myphr)
  PHR.df$Level = 'PHR'
  colnames(PHR.df)[3] = 'Level_Name' 
  
  ### METRO ###
  metro.temp<-cleaned.county %>%
    group_by(Metro_Area, Date) %>%
    mutate(metro_Cases_Daily = sum(Cases_Daily, na.rm = TRUE)) %>%
    mutate(metro_Cases_Cumulative = sum(Cases_Cumulative, na.rm=TRUE)) %>%
    mutate(metro_pop_DSHS = sum(Population_DSHS, na.rm=TRUE)) %>%
    dplyr::select(Date, Metro_Area, metro_Cases_Daily, metro_Cases_Cumulative, metro_pop_DSHS) %>%
    distinct()
    
  metro.df <- data.frame(Date = metro.temp$Date,
                         Level = 'metro',
                         Level_Name = metro.temp$Metro_Area,
                         Cases_Daily = metro.temp$metro_Cases_Daily,
                         Cases_Cumulative = metro.temp$metro_Cases_Cumulative,
                         Population_DSHS = metro.temp$metro_pop_DSHS)
  
  # drop NA dates
  metro.df <- subset(metro.df, !is.na(Date) & Date>=as.Date('2020-03-04'))
  
  ### TSA ###
  TSA.df <- ts.data.clean(mytsa)
  TSA.df <- subset(TSA.df, !is.na(Date) & Date >= as.Date('2020-03-04'))
  
  TSA.df$Level <- 'TSA'
  colnames(TSA.df)[2] <- 'Level_Name'
  
  
  ### STATE ###
  state.temp<-TSA.df %>%
    group_by(Date) %>%
    mutate(state_Cases_Daily = sum(Cases_Daily, na.rm = TRUE)) %>%
    mutate(state_Cases_Cumulative = sum(Cases_Cumulative, na.rm = TRUE)) %>%
    mutate(state_pop_DSHS = sum(Population_DSHS, na.rm = TRUE)) %>%
    mutate(state_hosp = sum(Hospitalizations_Total, na.rm=TRUE))%>%
    dplyr::select(Date, state_Cases_Daily, state_Cases_Cumulative, state_pop_DSHS, state_hosp) %>%
    distinct()
    
  state.df <- data.frame(Date=state.temp$Date,
                         Level = 'State',
                         Level_Name = 'Texas',
                         Cases_Daily=state.temp$state_Cases_Daily,
                         Cases_Cumulative=state.temp$state_Cases_Cumulative, 
                         Population_DSHS=state.temp$state_pop_DSHS,
                         Hospitalizations_Total=state.temp$state_hosp)
  
  state.df <- subset(state.df, Date>=as.Date('2020-03-04'))
  
  
  ### OUTPUT ###
  ts.df.out <- list(county=county.df, TSA=TSA.df, state=state.df, metro=metro.df, phr = PHR.df)
  return(ts.df.out)
}
# Compute forecast (UPDATE PREDICTION PERIOD [days] AS NEEDED)
covid.arima.forecast<-function(mydata, prediction.period = 10, mindate)
{
  maxdate <- max(mydata$Date)
  # mindate <- as.Date('2020-05-01')
  pred_start_label = format(mindate, format = '%m_%d')
  
  mydata = subset(mydata, Date >= mindate)
  model.length <- as.numeric(length(mydata$Date) + prediction.period)

  if(max(mydata$Cases_Daily>=100, na.rm = TRUE))
  {
    # arima requires cases to be a timeseries vector
    #convert daily cases to time series
    my.timeseries<-ts(mydata$Cases_Daily)
    print(mydata[1,2])
    #load package(pracma)
    library(pracma)
    my.timeseries<-movavg(my.timeseries,7,"s")
    #d=0 restricts first differencing to 0 so that daily cases aren't differenced
    
    arima.fit <- forecast::auto.arima(my.timeseries)
  
    # obtain diagnostic plots for ideal arima (p,d,q) selection 
    acf <- ggAcf(my.timeseries, lag.max=30)
    pacf <- ggPacf(my.timeseries, lag.max=30)
    ts.diagnostics <- grid.arrange(acf, pacf, nrow=2)
    ggsave(paste0('statistical-output/time-series/diagnostics/',
                  mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'),
           plot = ts.diagnostics)

    
    
    # save parameters from arima autofit
    p<- arima.fit$arma[1]          # autoregressive order 
    q<- arima.fit$arma[2]          # moving average order 
    d<-arima.fit$arma[6]           # differencing order from model
  
    # 10 day forecast, CI for lower and upper has confidence level 70% set by level =c(70,70)
    arima.forecast <- forecast::forecast(arima.fit, h = prediction.period, level=c(70,70))

    #return a dataframe of the arima model(Daily cases by date)
    arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
                            Cases_Raw = c(mydata$Cases_Daily, rep(NA, times = prediction.period)),
                            Cases_Daily = c(my.timeseries, arima.forecast[['mean']]),
                            CI_Lower = c(rep(NA, times = length(my.timeseries)),
                                         arima.forecast[['lower']][, 2]),
                            CI_Upper = c(rep(NA, times = length(my.timeseries)),
                                         arima.forecast[['upper']][, 2]))
                            # Order_AutoReg = c(rep(p, times = model.length)),
                            # Order_Moving_Avg = c(rep(q, times = model.length)),
                            # Differencing = c(rep(d, times = model.length)))
    
      # save prediction plot for preliminary review
      arima.plot <- ggplot(arima.out, aes(x=Date, y = Cases_Daily))+
                    geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), fill = "red", alpha = 0.5, size = 0.1)+
                    geom_line(color = "black", size = 1)+
                    labs(y = 'Daily Cases (7-Day Moving Average)', x = 'Date',
                         title = mydata$Level_Name[1]) +
                    scale_x_date(limits = c(mindate, maxdate + prediction.period),
                                 date_labels = '%b-%d', date_breaks = '1 week') + 
                    ggpubr::theme_pubr() +
                    theme(axis.text.x = element_text(angle = -90))
      
      ggsave(plot = arima.plot, paste0('statistical-output/time-series/plots/',
                                       mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'))    
    
    } else {
    # insufficient data catch: return NA values for predictions 
    arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
                            Cases_Raw = c(mydata$Cases_Daily, rep(NA, times = prediction.period)),
                            Cases_Daily = rep(NA, times = model.length),
                            CI_Lower = rep(NA, times = model.length),
                            CI_Upper =  rep(NA, times = model.length))
                            # Order_AutoReg = c(rep(NA, times = model.length)),
                            # Order_Moving_Avg = c(rep(NA, times = model.length)),
                            # Differencing = c(rep(NA, times= model.length)))
    }
  
  #replace CI lower limit with 0 if negative
  arima.out$CI_Lower <- ifelse(arima.out$CI_Lower>=0,arima.out$CI_Lower, 0)
  
  
  return(arima.out)
}

County

# apply arima across all counties
county.arima.output.1 <- nlme::gapply(county.Daily,
                                    FUN = covid.arima.forecast,
                                    groups = county.Daily$Level_Name,
                                    mindate = as.Date('2020-03-04'))
## [1] Anderson
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## 
## Attaching package: 'pracma'
## The following object is masked from 'package:mgcv':
## 
##     magic
## The following object is masked from 'package:purrr':
## 
##     cross
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Angelina
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Bell
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Bexar
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Brazoria
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Brazos
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Cameron
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Chambers
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Collin
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Dallas
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Denton
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Ector
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] El Paso
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Ellis
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Fort Bend
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Galveston
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Gregg
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Grimes
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Harris
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Hays
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Hidalgo
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Jefferson
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Johnson
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Jones
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Kaufman
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Lubbock
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Maverick
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] McLennan
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Medina
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Midland
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Montgomery
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Moore
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Nueces
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Orange
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Potter
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Randall
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Rusk
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Scurry
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Smith
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Starr
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Tarrant
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Titus
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Tom Green
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Travis
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Val Verde
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Victoria
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Walker
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Webb
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Williamson
## 254 Levels: Anderson Andrews Angelina Aransas Archer Armstrong ... Zavala
## Saving 7 x 5 in image
## Saving 7 x 5 in image

# bind list of dataframes to one dataframe
county.arima.df.1 <- data.table::rbindlist(county.arima.output.1, idcol = TRUE)
colnames(county.arima.df.1)[1] <- 'County'

# county.arima.output.2 <- nlme::gapply(county.Daily,
#                                     FUN = covid.arima.forecast,
#                                     groups = county.Daily$Level_Name,
#                                     mindate = as.Date('2020-06-03'))
# 
# 
# # bind list of dataframes to one dataframe
# county.arima.df.2 <- data.table::rbindlist(county.arima.output.2, idcol = TRUE)
# colnames(county.arima.df.2)[1] <- 'County'

Metro

# apply arima across both metro values
metro.arima.output.1 <- nlme::gapply(metro.Daily,
                                   FUN = covid.arima.forecast,
                                   groups = metro.Daily$Level_Name,
                                   mindate = as.Date('2020-03-04'))
## [1] metro
## Levels: metro
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] metro
## Levels: metro
## Saving 7 x 5 in image
## Saving 7 x 5 in image

# bind list of dataframes to one dataframe
metro.arima.df.1 <- data.table::rbindlist(metro.arima.output.1, idcol = TRUE)
colnames(metro.arima.df.1)[1] <- 'Metro_Area'
# 
# metro.arima.output.2 <- nlme::gapply(metro.Daily,
#                                    FUN = covid.arima.forecast,
#                                    groups = metro.Daily$Level_Name,
#                                    mindate = as.Date('2020-06-03'))
# 
# # bind list of dataframes to one dataframe
# metro.arima.df.2 <- data.table::rbindlist(metro.arima.output.2, idcol = TRUE)
# colnames(metro.arima.df.2)[1] <- 'Metro_Area'

TSA

# All data
TSA.arima.output.1 <- nlme::gapply(TSA.Daily,
                                   FUN = covid.arima.forecast,
                                   groups = TSA.Daily$Level_Name,
                                   mindate = as.Date('2020-03-04'))
## [1] A
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] B
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] D
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] E
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] F
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] G
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] H
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] I
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] J
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] K
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] L
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] M
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] N
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] O
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] P
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] Q
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] R
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] S
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] T
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] U
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] V
## Levels: A B C D E F G H I J K L M N O P Q R S T U V
## Saving 7 x 5 in image
## Saving 7 x 5 in image

TSA.arima.df.1 <- data.table::rbindlist(TSA.arima.output.1, idcol=TRUE)

colnames(TSA.arima.df.1)[1] <- 'TSA'
TSA.arima.df.1 <- merge(TSA.arima.df.1, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = c('TSA'))

# combine TSA 
TSA.arima.df.1$TSA <- paste0(TSA.arima.df.1$TSA, ' - ', TSA.arima.df.1$TSA_Name)
TSA.arima.df.1$TSA_Name <- NULL

PHR

PHR.arima.output <- nlme::gapply(PHR.Daily,
                                 FUN = covid.arima.forecast,
                                 groups = PHR.Daily$Level_Name,
                                 mindate = as.Date('2020-03-04'))
## [1] 2/3
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] 9/10
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] 11
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] 6/5S
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] 1
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] 8
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] 7
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10
## Saving 7 x 5 in image
## Saving 7 x 5 in image

## [1] 4/5N
## Levels: 1 11 2/3 4/5N 6/5S 7 8 9/10
## Saving 7 x 5 in image
## Saving 7 x 5 in image

PHR.arima.df <- data.table::rbindlist(PHR.arima.output, idcol=TRUE)
colnames(PHR.arima.df)[1] <- 'Level_Name'

PHR.arima.df <- merge(PHR.arima.df, unique(PHR.Daily[, c('PHR', 'Level_Name')]), by = c('Level_Name'))

# combine PHR 
PHR.arima.df$Level_Name <- paste0(PHR.arima.df$PHR, ' - ', PHR.arima.df$Level_Name)
PHR.arima.df$PHR <- NULL
colnames(PHR.arima.df)[1] <- 'PHR'

State

texas.arima.df.1 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-03-04'))
## Saving 7 x 5 in image
## Saving 7 x 5 in image

# texas.arima.df.2 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-06-03'))
#Create plot function, forecast.data is the dataframe containing true data
# and the forecast data. days.forecast is the number of days that are forecast
# in the dataframe (here we are forecasting 10 days, so it should be 10 unless we
# change this at some point down the road)
forecast.plot<-function(forecast.data, plot.title, y.label){
    mindate<-as.POSIXct("2020-03-04")
    maxdate<-as.POSIXct(max(forecast.data$Date))
  ts.plot<-ggplot(aes(x=as.POSIXct(Date), y=Cases_Daily), data=forecast.data)+
    geom_ribbon(aes(ymin=CI_Lower, ymax=CI_Upper), fill="grey50", size=0.1)+
    geom_line(color="blue", size=1)+
    geom_line(aes(x=as.POSIXct(Date), y=CI_Lower), color="grey", size=0.1)+
    geom_line(aes(x=as.POSIXct(Date),y=CI_Upper), color="grey", size=0.1)+
    scale_x_datetime(limits = c(mindate, maxdate))+
    xlab("Date")+
    ylab(y.label)+
    #Can use the following title if we are running using nlme for data frame
    #ggtitle(paste(forecast.data$.id,": TS Daily Cases + 10 Day Forcast",sep=""))
    ggtitle(plot.title)
ts.plot
}
######### View Output Table & Graph for TSA Q ##########
library(kableExtra)
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
#subset TSA Q arima data
TSA_Q.arima.df<-TSA.arima.df.1%>%subset(TSA.arima.df.1$TSA== "Q - Houston")
#apply the plot function to TSA Q arima data frame
TSA_Q.arima.df%>%kable(caption="Greater Houston ARIMA - Daily Cases")%>%kable_styling(full_width=FALSE)
Greater Houston ARIMA - Daily Cases
TSA Date Cases_Raw Cases_Daily CI_Lower CI_Upper
Q - Houston 2020-03-04 0 0.000000 NA NA
Q - Houston 2020-03-05 0 0.000000 NA NA
Q - Houston 2020-03-06 5 1.666667 NA NA
Q - Houston 2020-03-07 NA NA NA NA
Q - Houston 2020-03-08 NA NA NA NA
Q - Houston 2020-03-09 6 NA NA NA
Q - Houston 2020-03-10 0 NA NA NA
Q - Houston 2020-03-11 0 NA NA NA
Q - Houston 2020-03-12 3 NA NA NA
Q - Houston 2020-03-13 0 NA NA NA
Q - Houston 2020-03-14 NA NA NA NA
Q - Houston 2020-03-15 9 NA NA NA
Q - Houston 2020-03-16 0 NA NA NA
Q - Houston 2020-03-17 0 NA NA NA
Q - Houston 2020-03-18 0 NA NA NA
Q - Houston 2020-03-19 2 NA NA NA
Q - Houston 2020-03-20 17 NA NA NA
Q - Houston 2020-03-21 2 4.285714 NA NA
Q - Houston 2020-03-22 1 3.142857 NA NA
Q - Houston 2020-03-23 0 3.142857 NA NA
Q - Houston 2020-03-24 60 11.714286 NA NA
Q - Houston 2020-03-25 109 27.285714 NA NA
Q - Houston 2020-03-26 72 37.285714 NA NA
Q - Houston 2020-03-27 60 43.428571 NA NA
Q - Houston 2020-03-28 48 50.000000 NA NA
Q - Houston 2020-03-29 242 84.428571 NA NA
Q - Houston 2020-03-30 102 99.000000 NA NA
Q - Houston 2020-03-31 52 97.857143 NA NA
Q - Houston 2020-04-01 190 109.428571 NA NA
Q - Houston 2020-04-02 227 131.571429 NA NA
Q - Houston 2020-04-03 138 142.714286 NA NA
Q - Houston 2020-04-04 212 166.142857 NA NA
Q - Houston 2020-04-05 204 160.714286 NA NA
Q - Houston 2020-04-06 138 165.857143 NA NA
Q - Houston 2020-04-07 473 226.000000 NA NA
Q - Houston 2020-04-08 477 267.000000 NA NA
Q - Houston 2020-04-09 231 267.571429 NA NA
Q - Houston 2020-04-10 781 359.428571 NA NA
Q - Houston 2020-04-11 275 368.428571 NA NA
Q - Houston 2020-04-12 333 386.857143 NA NA
Q - Houston 2020-04-13 76 378.000000 NA NA
Q - Houston 2020-04-14 163 333.714286 NA NA
Q - Houston 2020-04-15 292 307.285714 NA NA
Q - Houston 2020-04-16 236 308.000000 NA NA
Q - Houston 2020-04-17 267 234.571429 NA NA
Q - Houston 2020-04-18 277 234.857143 NA NA
Q - Houston 2020-04-19 211 217.428571 NA NA
Q - Houston 2020-04-20 194 234.285714 NA NA
Q - Houston 2020-04-21 216 241.857143 NA NA
Q - Houston 2020-04-22 201 228.857143 NA NA
Q - Houston 2020-04-23 195 223.000000 NA NA
Q - Houston 2020-04-24 181 210.714286 NA NA
Q - Houston 2020-04-25 214 201.714286 NA NA
Q - Houston 2020-04-26 171 196.000000 NA NA
Q - Houston 2020-04-27 161 191.285714 NA NA
Q - Houston 2020-04-28 139 180.285714 NA NA
Q - Houston 2020-04-29 219 182.857143 NA NA
Q - Houston 2020-04-30 276 194.428571 NA NA
Q - Houston 2020-05-01 275 207.857143 NA NA
Q - Houston 2020-05-02 251 213.142857 NA NA
Q - Houston 2020-05-03 243 223.428571 NA NA
Q - Houston 2020-05-04 182 226.428571 NA NA
Q - Houston 2020-05-05 136 226.000000 NA NA
Q - Houston 2020-05-06 257 231.428571 NA NA
Q - Houston 2020-05-07 181 217.857143 NA NA
Q - Houston 2020-05-08 213 209.000000 NA NA
Q - Houston 2020-05-09 263 210.714286 NA NA
Q - Houston 2020-05-10 224 208.000000 NA NA
Q - Houston 2020-05-11 92 195.142857 NA NA
Q - Houston 2020-05-12 344 224.857143 NA NA
Q - Houston 2020-05-13 310 232.428571 NA NA
Q - Houston 2020-05-14 240 240.857143 NA NA
Q - Houston 2020-05-15 243 245.142857 NA NA
Q - Houston 2020-05-16 332 255.000000 NA NA
Q - Houston 2020-05-17 141 243.142857 NA NA
Q - Houston 2020-05-18 343 279.000000 NA NA
Q - Houston 2020-05-19 214 260.428571 NA NA
Q - Houston 2020-05-20 391 272.000000 NA NA
Q - Houston 2020-05-21 167 261.571429 NA NA
Q - Houston 2020-05-22 243 261.571429 NA NA
Q - Houston 2020-05-23 288 255.285714 NA NA
Q - Houston 2020-05-24 275 274.428571 NA NA
Q - Houston 2020-05-25 171 249.857143 NA NA
Q - Houston 2020-05-26 126 237.285714 NA NA
Q - Houston 2020-05-27 561 261.571429 NA NA
Q - Houston 2020-05-28 463 303.857143 NA NA
Q - Houston 2020-05-29 262 306.571429 NA NA
Q - Houston 2020-05-30 379 319.571429 NA NA
Q - Houston 2020-05-31 752 387.714286 NA NA
Q - Houston 2020-06-01 83 375.142857 NA NA
Q - Houston 2020-06-02 526 432.285714 NA NA
Q - Houston 2020-06-03 603 438.285714 NA NA
Q - Houston 2020-06-04 375 425.714286 NA NA
Q - Houston 2020-06-05 484 457.428571 NA NA
Q - Houston 2020-06-06 446 467.000000 NA NA
Q - Houston 2020-06-07 489 429.428571 NA NA
Q - Houston 2020-06-08 163 440.857143 NA NA
Q - Houston 2020-06-09 416 425.142857 NA NA
Q - Houston 2020-06-10 424 399.571429 NA NA
Q - Houston 2020-06-11 455 411.000000 NA NA
Q - Houston 2020-06-12 385 396.857143 NA NA
Q - Houston 2020-06-13 427 394.142857 NA NA
Q - Houston 2020-06-14 415 383.571429 NA NA
Q - Houston 2020-06-15 221 391.857143 NA NA
Q - Houston 2020-06-16 572 414.142857 NA NA
Q - Houston 2020-06-17 606 440.142857 NA NA
Q - Houston 2020-06-18 708 476.285714 NA NA
Q - Houston 2020-06-19 544 499.000000 NA NA
Q - Houston 2020-06-20 1443 644.142857 NA NA
Q - Houston 2020-06-21 1224 759.714286 NA NA
Q - Houston 2020-06-22 305 771.714286 NA NA
Q - Houston 2020-06-23 2187 1002.428571 NA NA
Q - Houston 2020-06-24 1581 1141.714286 NA NA
Q - Houston 2020-06-25 1575 1265.571429 NA NA
Q - Houston 2020-06-26 1438 1393.285714 NA NA
Q - Houston 2020-06-27 1604 1416.285714 NA NA
Q - Houston 2020-06-28 1019 1387.000000 NA NA
Q - Houston 2020-06-29 142 1363.714286 NA NA
Q - Houston 2020-06-30 1569 1275.428571 NA NA
Q - Houston 2020-07-01 953 1185.714286 NA NA
Q - Houston 2020-07-02 1713 1205.428571 NA NA
Q - Houston 2020-07-03 1554 1222.000000 NA NA
Q - Houston 2020-07-04 1356 1186.571429 NA NA
Q - Houston 2020-07-05 594 1125.857143 NA NA
Q - Houston 2020-07-06 734 1210.428571 NA NA
Q - Houston 2020-07-07 1580 1212.000000 NA NA
Q - Houston 2020-07-08 1785 1330.857143 NA NA
Q - Houston 2020-07-09 996 1228.428571 NA NA
Q - Houston 2020-07-10 1250 1185.000000 NA NA
Q - Houston 2020-07-11 1383 1188.857143 NA NA
Q - Houston 2020-07-12 2137 1409.285714 NA NA
Q - Houston 2020-07-13 1524 1522.142857 NA NA
Q - Houston 2020-07-14 2595 1667.142857 NA NA
Q - Houston 2020-07-15 2253 1734.000000 NA NA
Q - Houston 2020-07-16 2314 1922.285714 NA NA
Q - Houston 2020-07-17 1935 2020.142857 NA NA
Q - Houston 2020-07-18 NA 2117.348488 2068.476 2166.221
Q - Houston 2020-07-19 NA 2189.896663 2108.721 2271.072
Q - Houston 2020-07-20 NA 2254.952757 2137.989 2371.916
Q - Houston 2020-07-21 NA 2311.115242 2159.979 2462.251
Q - Houston 2020-07-22 NA 2362.612133 2178.114 2547.111
Q - Houston 2020-07-23 NA 2410.318713 2193.772 2626.866
Q - Houston 2020-07-24 NA 2455.638505 2208.168 2703.109
Q - Houston 2020-07-25 NA 2499.224726 2221.898 2776.551
Q - Houston 2020-07-26 NA 2541.650889 2235.388 2847.914
Q - Houston 2020-07-27 NA 2583.263871 2248.863 2917.665
#nlme::gapply(TSA.arima.df, FUN = forecast.plot, groups=TSA.daily$Level_Name)
forecast.plot(TSA_Q.arima.df, "TSA Q - Greater Houston Daily Cases", "Daily Cases")
## Warning: Removed 1 rows containing missing values (geom_path).
## Warning: Removed 136 rows containing missing values (geom_path).

## Warning: Removed 136 rows containing missing values (geom_path).

######### View Output Table & Graph for Texas ##########
texas.arima.df.1%>%kable(caption="Texas Arima - Daily Cases")%>%kable_styling(full_width=FALSE)
Texas Arima - Daily Cases
Date Cases_Raw Cases_Daily CI_Lower CI_Upper
2020-03-04 0 0.000000 NA NA
2020-03-05 0 0.000000 NA NA
2020-03-06 5 1.666667 NA NA
2020-03-07 0 1.250000 NA NA
2020-03-08 0 1.000000 NA NA
2020-03-09 7 2.000000 NA NA
2020-03-10 3 2.142857 NA NA
2020-03-11 3 2.571429 NA NA
2020-03-12 4 3.142857 NA NA
2020-03-13 0 2.428571 NA NA
2020-03-14 0 2.428571 NA NA
2020-03-15 34 7.285714 NA NA
2020-03-16 0 6.285714 NA NA
2020-03-17 7 6.857143 NA NA
2020-03-18 19 9.142857 NA NA
2020-03-19 26 12.285714 NA NA
2020-03-20 67 21.857143 NA NA
2020-03-21 60 30.428571 NA NA
2020-03-22 28 29.571429 NA NA
2020-03-23 24 33.000000 NA NA
2020-03-24 425 92.714286 NA NA
2020-03-25 263 127.571429 NA NA
2020-03-26 419 183.714286 NA NA
2020-03-27 337 222.285714 NA NA
2020-03-28 317 259.000000 NA NA
2020-03-29 504 327.000000 NA NA
2020-03-30 322 369.571429 NA NA
2020-03-31 392 364.857143 NA NA
2020-04-01 730 431.571429 NA NA
2020-04-02 669 467.285714 NA NA
2020-04-03 659 513.285714 NA NA
2020-04-04 788 580.571429 NA NA
2020-04-05 681 605.857143 NA NA
2020-04-06 484 629.000000 NA NA
2020-04-07 988 714.142857 NA NA
2020-04-08 1092 765.857143 NA NA
2020-04-09 877 795.571429 NA NA
2020-04-10 1441 907.285714 NA NA
2020-04-11 890 921.857143 NA NA
2020-04-12 923 956.428571 NA NA
2020-04-13 422 947.571429 NA NA
2020-04-14 718 909.000000 NA NA
2020-04-15 868 877.000000 NA NA
2020-04-16 963 889.285714 NA NA
2020-04-17 916 814.285714 NA NA
2020-04-18 889 814.142857 NA NA
2020-04-19 663 777.000000 NA NA
2020-04-20 535 793.142857 NA NA
2020-04-21 738 796.000000 NA NA
2020-04-22 873 796.714286 NA NA
2020-04-23 876 784.285714 NA NA
2020-04-24 862 776.571429 NA NA
2020-04-25 968 787.857143 NA NA
2020-04-26 859 815.857143 NA NA
2020-04-27 666 834.571429 NA NA
2020-04-28 874 854.000000 NA NA
2020-04-29 885 855.714286 NA NA
2020-04-30 1033 878.142857 NA NA
2020-05-01 1142 918.142857 NA NA
2020-05-02 1293 964.571429 NA NA
2020-05-03 1026 988.428571 NA NA
2020-05-04 784 1005.285714 NA NA
2020-05-05 1037 1028.571429 NA NA
2020-05-06 1053 1052.571429 NA NA
2020-05-07 1135 1067.142857 NA NA
2020-05-08 1219 1078.142857 NA NA
2020-05-09 1251 1072.142857 NA NA
2020-05-10 1009 1069.714286 NA NA
2020-05-11 1000 1100.571429 NA NA
2020-05-12 1179 1120.857143 NA NA
2020-05-13 1355 1164.000000 NA NA
2020-05-14 1448 1208.714286 NA NA
2020-05-15 1347 1227.000000 NA NA
2020-05-16 1801 1305.571429 NA NA
2020-05-17 785 1273.571429 NA NA
2020-05-18 909 1260.571429 NA NA
2020-05-19 1219 1266.285714 NA NA
2020-05-20 1411 1274.285714 NA NA
2020-05-21 945 1202.428571 NA NA
2020-05-22 1241 1187.285714 NA NA
2020-05-23 1060 1081.428571 NA NA
2020-05-24 839 1089.142857 NA NA
2020-05-25 623 1048.285714 NA NA
2020-05-26 620 962.714286 NA NA
2020-05-27 1361 955.571429 NA NA
2020-05-28 1855 1085.571429 NA NA
2020-05-29 1230 1084.000000 NA NA
2020-05-30 1333 1123.000000 NA NA
2020-05-31 1949 1281.571429 NA NA
2020-06-01 593 1277.285714 NA NA
2020-06-02 1688 1429.857143 NA NA
2020-06-03 1703 1478.714286 NA NA
2020-06-04 1650 1449.428571 NA NA
2020-06-05 1693 1515.571429 NA NA
2020-06-06 1942 1602.571429 NA NA
2020-06-07 1426 1527.857143 NA NA
2020-06-08 638 1534.285714 NA NA
2020-06-09 1853 1557.857143 NA NA
2020-06-10 2509 1673.000000 NA NA
2020-06-11 1826 1698.142857 NA NA
2020-06-12 2097 1755.857143 NA NA
2020-06-13 2331 1811.428571 NA NA
2020-06-14 1843 1871.000000 NA NA
2020-06-15 1254 1959.000000 NA NA
2020-06-16 4098 2279.714286 NA NA
2020-06-17 3140 2369.857143 NA NA
2020-06-18 3516 2611.285714 NA NA
2020-06-19 3454 2805.142857 NA NA
2020-06-20 4430 3105.000000 NA NA
2020-06-21 3866 3394.000000 NA NA
2020-06-22 3280 3683.428571 NA NA
2020-06-23 5522 3886.857143 NA NA
2020-06-24 5551 4231.285714 NA NA
2020-06-25 5996 4585.571429 NA NA
2020-06-26 5707 4907.428571 NA NA
2020-06-27 5742 5094.857143 NA NA
2020-06-28 5357 5307.857143 NA NA
2020-06-29 4288 5451.857143 NA NA
2020-06-30 6975 5659.428571 NA NA
2020-07-01 8076 6020.142857 NA NA
2020-07-02 7915 6294.285714 NA NA
2020-07-03 7555 6558.285714 NA NA
2020-07-04 8258 6917.714286 NA NA
2020-07-05 3449 6645.142857 NA NA
2020-07-06 5319 6792.428571 NA NA
2020-07-07 10028 7228.571429 NA NA
2020-07-08 9979 7500.428571 NA NA
2020-07-09 9782 7767.142857 NA NA
2020-07-10 9765 8082.857143 NA NA
2020-07-11 10351 8381.857143 NA NA
2020-07-12 8196 9060.000000 NA NA
2020-07-13 5655 9108.000000 NA NA
2020-07-14 10745 9210.428571 NA NA
2020-07-15 9550 9149.142857 NA NA
2020-07-16 10291 9221.857143 NA NA
2020-07-17 14916 9957.714286 NA NA
2020-07-18 NA 10228.153935 10117.32 10338.99
2020-07-19 NA 10498.593584 10328.31 10668.88
2020-07-20 NA 10769.033233 10543.61 10994.46
2020-07-21 NA 11039.472882 10759.42 11319.52
2020-07-22 NA 11309.912531 10974.49 11645.33
2020-07-23 NA 11580.352181 11188.26 11972.44
2020-07-24 NA 11850.791830 11400.47 12301.11
2020-07-25 NA 12121.231479 11610.98 12631.48
2020-07-26 NA 12391.671128 11819.73 12963.61
2020-07-27 NA 12662.110777 12026.71 13297.51
forecast.plot(texas.arima.df.1, "Texas Daily Cases", "Daily Cases")
## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 136 rows containing missing values (geom_path).

## Warning: Removed 136 rows containing missing values (geom_path).

export

write.csv(texas.arima.df.1, 'statistical-output/time-series/state_case_timeseries.csv', row.names = FALSE)
write.csv(TSA.arima.df.1, 'statistical-output/time-series/tsa_case_timeseries.csv', row.names = FALSE)
write.csv(county.arima.df.1, 'statistical-output/time-series/county_case_timeseries.csv', row.names = FALSE)
write.csv(metro.arima.df.1, 'statistical-output/time-series/metro_case_timeseries.csv', row.names = FALSE)
write.csv(PHR.arima.df, 'statistical-output/time-series/phr_case_timeseries.csv', row.names = FALSE)

# write.csv(texas.arima.df.2, 'statistical-output/time-series/state_timeseries_06_03.csv', row.names = FALSE)
# write.csv(TSA.arima.df.2, 'statistical-output/time-series/tsa_timeseries_06_03.csv', row.names = FALSE)
# write.csv(county.arima.df.2, 'statistical-output/time-series/county_timeseries_06_03.csv', row.names = FALSE)
# write.csv(metro.arima.df.2, 'statistical-output/time-series/metro_timeseries_06_03.csv', row.names = FALSE)

grouping

colnames(county.arima.df.1)[1] <- 'Level'
colnames(metro.arima.df.1)[1] <- 'Level'
colnames(TSA.arima.df.1)[1] <- 'Level'
texas.arima.df.1$Level <- 'Texas'
colnames(PHR.arima.df)[1] <- 'Level'

county.arima.df.1$Level_Type = 'County'
metro.arima.df.1$Level_Type = 'Metro'
TSA.arima.df.1$Level_Type = 'TSA'
PHR.arima.df$Level_Type = 'PHR'
texas.arima.df.1$Level_Type = 'State'


combined.arima.df.1 <- rbind(county.arima.df.1, metro.arima.df.1, TSA.arima.df.1, texas.arima.df.1, PHR.arima.df)
write.csv(combined.arima.df.1, 'statistical-output/time-series/stacked_case_timeseries.csv',
          row.names = FALSE)

write.csv(combined.arima.df.1, 'tableau/stacked_case_timeseries.csv',
          row.names = FALSE)
# 
# combined.arima.df.2 <- rbind(county.arima.df.2, metro.arima.df.2, TSA.arima.df.2, texas.arima.df.2)
# write.csv(combined.arima.df.2, 'statistical-output/time-series/stacked_timeseries_06_03.csv',
#           row.names = FALSE)
# 
# write.csv(combined.arima.df.2, 'tableau/stacked_timeseries_06_03.csv',
#           row.names = FALSE)

Hospitalization Time Series

# Compute forecast (UPDATE PREDICTION PERIOD [days] AS NEEDED)
covid.arima.forecast<-function(mydata, prediction.period = 10, mindate)
{
  maxdate <- max(mydata$Date)
  # mindate <- as.Date('2020-05-01')
  pred_start_label = format(mindate, format = '%m_%d')
  
  mydata = subset(mydata, Date >= mindate)
  model.length <- as.numeric(length(mydata$Date) + prediction.period)

  if(max(mydata$Hospitalizations_Total>=100, na.rm = TRUE))
  {

    
    # arima requires cases to be a timeseries vector
    #convert daily cases to time series
    my.timeseries<-ts(mydata$Hospitalizations_Total)
    #my.timeseries <- rollmeanr(my.timeseries, k=7, na.pad=TRUE, align = 'right')
    
    #load package(pracma)
    library(pracma)
    my.timeseries<-movavg(my.timeseries,7,"s")
    #d=0 restricts first differencing to 0 so that daily cases aren't differenced
    
    arima.fit <- forecast::auto.arima(my.timeseries)
  
    # obtain diagnostic plots for ideal arima (p,d,q) selection 
    acf <- ggAcf(my.timeseries, lag.max=30)
    pacf <- ggPacf(my.timeseries, lag.max=30)
    ts.diagnostics <- grid.arrange(acf, pacf, nrow=2)
    ggsave(paste0('statistical-output/time-series/diagnostics/',
                  mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'),
           plot = ts.diagnostics)

    
    
    # save parameters from arima autofit
    p<- arima.fit$arma[1]          # autoregressive order 
    q<- arima.fit$arma[2]          # moving average order 
    d<-arima.fit$arma[6]           # differencing order from model
  
    # 10 day forecast, CI for lower and upper has confidence level 70% set by level =c(70,70)
    arima.forecast <- forecast::forecast(arima.fit, h = prediction.period, level=c(70,70))

    #return a dataframe of the arima model(Daily cases by date)
    arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
                            Cases_Raw = c(mydata$Hospitalizations_Total, rep(NA, times = prediction.period)),
                            Hospitalizations_Total = c(my.timeseries, arima.forecast[['mean']]),
                            CI_Lower = c(rep(NA, times = length(my.timeseries)),
                                         arima.forecast[['lower']][, 2]),
                            CI_Upper = c(rep(NA, times = length(my.timeseries)),
                                         arima.forecast[['upper']][, 2]))
                            # Order_AutoReg = c(rep(p, times = model.length)),
                            # Order_Moving_Avg = c(rep(q, times = model.length)),
                            # Differencing = c(rep(d, times = model.length)))
    
      # save prediction plot for preliminary review
      arima.plot <- ggplot(arima.out, aes(x=Date, y = Hospitalizations_Total))+
                    geom_ribbon(aes(ymin = CI_Lower, ymax = CI_Upper), fill = "red", alpha = 0.5, size = 0.1)+
                    geom_line(color = "black", size = 1)+
                    labs(y = 'Daily Cases (7-Day Moving Average)', x = 'Date',
                         title = mydata$Level_Name[1]) +
                    scale_x_date(limits = c(mindate, maxdate + prediction.period),
                                 date_labels = '%b-%d', date_breaks = '1 week') + 
                    ggpubr::theme_pubr() +
                    theme(axis.text.x = element_text(angle = -90))
      
      ggsave(plot = arima.plot, paste0('statistical-output/time-series/plots/',
                                       mydata$Level[1],'/', mydata$Level_Name[1], pred_start_label, '.png'))    
    
    } else {
    # insufficient data catch: return NA values for predictions 
    arima.out <- data.frame(Date = seq(mindate, maxdate + prediction.period, by = 'days'),
                            Cases_Raw = c(mydata$Hospitalizations_Total, rep(NA, times = prediction.period)),
                            Hospitalizations_Total = rep(NA, times = model.length),
                            CI_Lower = rep(NA, times = model.length),
                            CI_Upper =  rep(NA, times = model.length))
                            # Order_AutoReg = c(rep(NA, times = model.length)),
                            # Order_Moving_Avg = c(rep(NA, times = model.length)),
                            # Differencing = c(rep(NA, times= model.length)))
    }
  
  #replace CI lower limit with 0 if negative
  arima.out$CI_Lower <- ifelse(arima.out$CI_Lower>=0,arima.out$CI_Lower, 0)
  
  
  return(arima.out)
}

TSA - Hospitalization

# All data
TSA.arima.output.1 <- nlme::gapply(TSA.Daily,
                                   FUN = covid.arima.forecast,
                                   groups = TSA.Daily$Level_Name,
                                   mindate = as.Date('2020-03-04'))
## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

## Saving 7 x 5 in image
## Saving 7 x 5 in image
## Warning: Removed 45 rows containing missing values (geom_path).

TSA.arima.df.1 <- data.table::rbindlist(TSA.arima.output.1, idcol=TRUE)

colnames(TSA.arima.df.1)[1] <- 'TSA'
TSA.arima.df.1 <- merge(TSA.arima.df.1, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = c('TSA'))

# combine TSA 
TSA.arima.df.1$TSA <- paste0(TSA.arima.df.1$TSA, ' - ', TSA.arima.df.1$TSA_Name)
TSA.arima.df.1$TSA_Name <- NULL

State - hospitalization

texas.arima.df.1 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-03-04'))
## Saving 7 x 5 in image
## Saving 7 x 5 in image

# texas.arima.df.2 <- covid.arima.forecast(state.Daily, mindate = as.Date('2020-06-03'))

Hospitalization Write to CSV

tsa.hosp.arima<-write.csv(TSA.arima.df.1, "statistical-output/time-series/tsa_hosp_timeseries.csv", row.names=FALSE)
texas.hosp.arima<-write.csv(texas.arima.df.1, "statistical-output/time-series/state_hosp_timeseries.csv", row.names=FALSE)

Stacking

TSA_hosp_out = TSA.arima.df.1
colnames(TSA_hosp_out)[1] = 'Level'
texas.arima.df.1$Level = 'Texas'


TSA_hosp_out$Level_Type = 'TSA'
texas.arima.df.1$Level_Type = 'State'

combined.arima.df.1 <- rbind(TSA_hosp_out, texas.arima.df.1)
write.csv(combined.arima.df.1, 'statistical-output/time-series/stacked_hosp_timeseries.csv',
          row.names = FALSE)

write.csv(combined.arima.df.1, 'tableau/stacked_hosp_timeseries.csv',
          row.names = FALSE)

Hospitalization Forecast Plots

#Create plot function, forecast.data is the dataframe containing true data
# and the forecast data. days.forecast is the number of days that are forecast
# in the dataframe (here we are forecasting 10 days, so it should be 10 unless we
# change this at some point down the road)
forecast.plot<-function(forecast.data, plot.title, y.label){
    mindate<-as.POSIXct("2020-03-04")
    maxdate<-as.POSIXct(max(forecast.data$Date))
  ts.plot<-ggplot(aes(x=as.POSIXct(Date), y=Hospitalizations_Total), data=forecast.data)+
    geom_ribbon(aes(ymin=CI_Lower, ymax=CI_Upper), fill="grey50", size=0.1)+
    geom_line(color="blue", size=1)+
    geom_line(aes(x=as.POSIXct(Date), y=CI_Lower), color="grey", size=0.1)+
    geom_line(aes(x=as.POSIXct(Date),y=CI_Upper), color="grey", size=0.1)+
    scale_x_datetime(limits = c(mindate, maxdate))+
    xlab("Date")+
    ylab(y.label)+
    #Can use the following title if we are running using nlme for data frame
    #ggtitle(paste(forecast.data$.id,": TS Daily Cases + 10 Day Forcast",sep=""))
    ggtitle(plot.title)
ts.plot
}
######### View Output Table & Graph for TSA Q ##########
library(kableExtra)
#subset TSA Q arima data
TSA_Q.arima.df<-TSA.arima.df.1%>%subset(TSA.arima.df.1$TSA == "Q - Houston")
#apply the plot function to TSA Q arima data frame
TSA_Q.arima.df%>%kable(caption="Greater Houston ARIMA - Daily Hosp")%>%kable_styling(full_width=FALSE)
Greater Houston ARIMA - Daily Hosp
TSA Date Cases_Raw Hospitalizations_Total CI_Lower CI_Upper
Q - Houston 2020-03-04 NA NA NA NA
Q - Houston 2020-03-05 NA NA NA NA
Q - Houston 2020-03-06 NA NA NA NA
Q - Houston 2020-03-07 NA NA NA NA
Q - Houston 2020-03-08 NA NA NA NA
Q - Houston 2020-03-09 NA NA NA NA
Q - Houston 2020-03-10 NA NA NA NA
Q - Houston 2020-03-11 NA NA NA NA
Q - Houston 2020-03-12 NA NA NA NA
Q - Houston 2020-03-13 NA NA NA NA
Q - Houston 2020-03-14 NA NA NA NA
Q - Houston 2020-03-15 NA NA NA NA
Q - Houston 2020-03-16 NA NA NA NA
Q - Houston 2020-03-17 NA NA NA NA
Q - Houston 2020-03-18 NA NA NA NA
Q - Houston 2020-03-19 NA NA NA NA
Q - Houston 2020-03-20 NA NA NA NA
Q - Houston 2020-03-21 NA NA NA NA
Q - Houston 2020-03-22 NA NA NA NA
Q - Houston 2020-03-23 NA NA NA NA
Q - Houston 2020-03-24 NA NA NA NA
Q - Houston 2020-03-25 NA NA NA NA
Q - Houston 2020-03-26 NA NA NA NA
Q - Houston 2020-03-27 NA NA NA NA
Q - Houston 2020-03-28 NA NA NA NA
Q - Houston 2020-03-29 NA NA NA NA
Q - Houston 2020-03-30 NA NA NA NA
Q - Houston 2020-03-31 NA NA NA NA
Q - Houston 2020-04-01 NA NA NA NA
Q - Houston 2020-04-02 NA NA NA NA
Q - Houston 2020-04-03 NA NA NA NA
Q - Houston 2020-04-04 NA NA NA NA
Q - Houston 2020-04-05 NA NA NA NA
Q - Houston 2020-04-06 NA NA NA NA
Q - Houston 2020-04-07 NA NA NA NA
Q - Houston 2020-04-08 NA NA NA NA
Q - Houston 2020-04-09 NA NA NA NA
Q - Houston 2020-04-10 NA NA NA NA
Q - Houston 2020-04-11 NA NA NA NA
Q - Houston 2020-04-12 15 NA NA NA
Q - Houston 2020-04-13 27 NA NA NA
Q - Houston 2020-04-14 37 NA NA NA
Q - Houston 2020-04-15 45 NA NA NA
Q - Houston 2020-04-16 52 NA NA NA
Q - Houston 2020-04-17 59 NA NA NA
Q - Houston 2020-04-18 64 42.71429 NA NA
Q - Houston 2020-04-19 72 50.85714 NA NA
Q - Houston 2020-04-20 15 49.14286 NA NA
Q - Houston 2020-04-21 15 46.00000 NA NA
Q - Houston 2020-04-22 86 51.85714 NA NA
Q - Houston 2020-04-23 91 57.42857 NA NA
Q - Houston 2020-04-24 97 62.85714 NA NA
Q - Houston 2020-04-25 102 68.28571 NA NA
Q - Houston 2020-04-26 107 73.28571 NA NA
Q - Houston 2020-04-27 110 86.85714 NA NA
Q - Houston 2020-04-28 114 101.00000 NA NA
Q - Houston 2020-04-29 117 105.42857 NA NA
Q - Houston 2020-04-30 122 109.85714 NA NA
Q - Houston 2020-05-01 126 114.00000 NA NA
Q - Houston 2020-05-02 130 118.00000 NA NA
Q - Houston 2020-05-03 134 121.85714 NA NA
Q - Houston 2020-05-04 137 125.71429 NA NA
Q - Houston 2020-05-05 141 129.57143 NA NA
Q - Houston 2020-05-06 117 129.57143 NA NA
Q - Houston 2020-05-07 149 133.42857 NA NA
Q - Houston 2020-05-08 152 137.14286 NA NA
Q - Houston 2020-05-09 154 140.57143 NA NA
Q - Houston 2020-05-10 154 143.42857 NA NA
Q - Houston 2020-05-11 161 146.85714 NA NA
Q - Houston 2020-05-12 117 143.42857 NA NA
Q - Houston 2020-05-13 70 136.71429 NA NA
Q - Houston 2020-05-14 166 139.14286 NA NA
Q - Houston 2020-05-15 168 141.42857 NA NA
Q - Houston 2020-05-16 173 144.14286 NA NA
Q - Houston 2020-05-17 175 147.14286 NA NA
Q - Houston 2020-05-18 178 149.57143 NA NA
Q - Houston 2020-05-19 181 158.71429 NA NA
Q - Houston 2020-05-20 78 159.85714 NA NA
Q - Houston 2020-05-21 187 162.85714 NA NA
Q - Houston 2020-05-22 191 166.14286 NA NA
Q - Houston 2020-05-23 194 169.14286 NA NA
Q - Houston 2020-05-24 196 172.14286 NA NA
Q - Houston 2020-05-25 199 175.14286 NA NA
Q - Houston 2020-05-26 168 173.28571 NA NA
Q - Houston 2020-05-27 173 186.85714 NA NA
Q - Houston 2020-05-28 70 170.14286 NA NA
Q - Houston 2020-05-29 196 170.85714 NA NA
Q - Houston 2020-05-30 209 173.00000 NA NA
Q - Houston 2020-05-31 211 175.14286 NA NA
Q - Houston 2020-06-01 213 177.14286 NA NA
Q - Houston 2020-06-02 215 183.85714 NA NA
Q - Houston 2020-06-03 91 172.14286 NA NA
Q - Houston 2020-06-04 218 193.28571 NA NA
Q - Houston 2020-06-05 220 196.71429 NA NA
Q - Houston 2020-06-06 222 198.57143 NA NA
Q - Houston 2020-06-07 225 200.57143 NA NA
Q - Houston 2020-06-08 226 202.42857 NA NA
Q - Houston 2020-06-09 228 204.28571 NA NA
Q - Houston 2020-06-10 231 224.28571 NA NA
Q - Houston 2020-06-11 233 226.42857 NA NA
Q - Houston 2020-06-12 236 228.71429 NA NA
Q - Houston 2020-06-13 240 231.28571 NA NA
Q - Houston 2020-06-14 243 233.85714 NA NA
Q - Houston 2020-06-15 249 237.14286 NA NA
Q - Houston 2020-06-16 253 240.71429 NA NA
Q - Houston 2020-06-17 257 244.42857 NA NA
Q - Houston 2020-06-18 263 248.71429 NA NA
Q - Houston 2020-06-19 264 252.71429 NA NA
Q - Houston 2020-06-20 271 257.14286 NA NA
Q - Houston 2020-06-21 276 261.85714 NA NA
Q - Houston 2020-06-22 281 266.42857 NA NA
Q - Houston 2020-06-23 288 271.42857 NA NA
Q - Houston 2020-06-24 293 276.57143 NA NA
Q - Houston 2020-06-25 300 281.85714 NA NA
Q - Houston 2020-06-26 307 288.00000 NA NA
Q - Houston 2020-06-27 314 294.14286 NA NA
Q - Houston 2020-06-28 320 300.42857 NA NA
Q - Houston 2020-06-29 327 307.00000 NA NA
Q - Houston 2020-06-30 333 313.42857 NA NA
Q - Houston 2020-07-01 340 320.14286 NA NA
Q - Houston 2020-07-02 348 327.00000 NA NA
Q - Houston 2020-07-03 355 333.85714 NA NA
Q - Houston 2020-07-04 365 341.14286 NA NA
Q - Houston 2020-07-05 372 348.57143 NA NA
Q - Houston 2020-07-06 379 356.00000 NA NA
Q - Houston 2020-07-07 388 363.85714 NA NA
Q - Houston 2020-07-08 394 371.57143 NA NA
Q - Houston 2020-07-09 403 379.42857 NA NA
Q - Houston 2020-07-10 412 387.57143 NA NA
Q - Houston 2020-07-11 421 395.57143 NA NA
Q - Houston 2020-07-12 429 403.71429 NA NA
Q - Houston 2020-07-13 421 409.71429 NA NA
Q - Houston 2020-07-14 447 418.14286 NA NA
Q - Houston 2020-07-15 456 427.00000 NA NA
Q - Houston 2020-07-16 460 435.14286 NA NA
Q - Houston 2020-07-17 467 443.00000 NA NA
Q - Houston 2020-07-18 NA 450.63939 445.6334 455.6454
Q - Houston 2020-07-19 NA 458.32642 451.5160 465.1368
Q - Houston 2020-07-20 NA 466.00302 457.3059 474.7002
Q - Houston 2020-07-21 NA 473.68191 463.1484 484.2154
Q - Houston 2020-07-22 NA 481.36030 468.9656 493.7550
Q - Houston 2020-07-23 NA 489.03879 474.7496 503.3280
Q - Houston 2020-07-24 NA 496.71726 480.4901 512.9445
Q - Houston 2020-07-25 NA 504.39574 486.1827 522.6088
Q - Houston 2020-07-26 NA 512.07422 491.8249 532.3235
Q - Houston 2020-07-27 NA 519.75269 497.4155 542.0899
#nlme::gapply(TSA.arima.df, FUN = forecast.plot, groups=TSA.daily$Level_Name)
forecast.plot(TSA_Q.arima.df, "TSA Q - Greater Houston Daily Hosp", "Daily Hosp")
## Warning: Removed 45 rows containing missing values (geom_path).
## Warning: Removed 136 rows containing missing values (geom_path).

## Warning: Removed 136 rows containing missing values (geom_path).

######### View Output Table & Graph for Texas ##########
texas.arima.df.1%>%kable(caption="Texas Arima - Daily Hosp")%>%kable_styling(full_width=FALSE)
Texas Arima - Daily Hosp
Date Cases_Raw Hospitalizations_Total CI_Lower CI_Upper Level Level_Type
2020-03-04 0 0.00000 NA NA Texas State
2020-03-05 0 0.00000 NA NA Texas State
2020-03-06 0 0.00000 NA NA Texas State
2020-03-07 0 0.00000 NA NA Texas State
2020-03-08 0 0.00000 NA NA Texas State
2020-03-09 0 0.00000 NA NA Texas State
2020-03-10 0 0.00000 NA NA Texas State
2020-03-11 0 0.00000 NA NA Texas State
2020-03-12 0 0.00000 NA NA Texas State
2020-03-13 0 0.00000 NA NA Texas State
2020-03-14 0 0.00000 NA NA Texas State
2020-03-15 0 0.00000 NA NA Texas State
2020-03-16 0 0.00000 NA NA Texas State
2020-03-17 0 0.00000 NA NA Texas State
2020-03-18 0 0.00000 NA NA Texas State
2020-03-19 0 0.00000 NA NA Texas State
2020-03-20 0 0.00000 NA NA Texas State
2020-03-21 0 0.00000 NA NA Texas State
2020-03-22 0 0.00000 NA NA Texas State
2020-03-23 0 0.00000 NA NA Texas State
2020-03-24 0 0.00000 NA NA Texas State
2020-03-25 0 0.00000 NA NA Texas State
2020-03-26 0 0.00000 NA NA Texas State
2020-03-27 0 0.00000 NA NA Texas State
2020-03-28 0 0.00000 NA NA Texas State
2020-03-29 0 0.00000 NA NA Texas State
2020-03-30 0 0.00000 NA NA Texas State
2020-03-31 0 0.00000 NA NA Texas State
2020-04-01 0 0.00000 NA NA Texas State
2020-04-02 0 0.00000 NA NA Texas State
2020-04-03 0 0.00000 NA NA Texas State
2020-04-04 0 0.00000 NA NA Texas State
2020-04-05 0 0.00000 NA NA Texas State
2020-04-06 0 0.00000 NA NA Texas State
2020-04-07 0 0.00000 NA NA Texas State
2020-04-08 0 0.00000 NA NA Texas State
2020-04-09 0 0.00000 NA NA Texas State
2020-04-10 0 0.00000 NA NA Texas State
2020-04-11 0 0.00000 NA NA Texas State
2020-04-12 229 32.71429 NA NA Texas State
2020-04-13 439 95.42857 NA NA Texas State
2020-04-14 547 173.57143 NA NA Texas State
2020-04-15 581 256.57143 NA NA Texas State
2020-04-16 570 338.00000 NA NA Texas State
2020-04-17 673 434.14286 NA NA Texas State
2020-04-18 827 552.28571 NA NA Texas State
2020-04-19 704 620.14286 NA NA Texas State
2020-04-20 739 663.00000 NA NA Texas State
2020-04-21 916 715.71429 NA NA Texas State
2020-04-22 1117 792.28571 NA NA Texas State
2020-04-23 1129 872.14286 NA NA Texas State
2020-04-24 1092 932.00000 NA NA Texas State
2020-04-25 1110 972.42857 NA NA Texas State
2020-04-26 1076 1025.57143 NA NA Texas State
2020-04-27 1172 1087.42857 NA NA Texas State
2020-04-28 1129 1117.85714 NA NA Texas State
2020-04-29 1278 1140.85714 NA NA Texas State
2020-04-30 1150 1143.85714 NA NA Texas State
2020-05-01 1283 1171.14286 NA NA Texas State
2020-05-02 1091 1168.42857 NA NA Texas State
2020-05-03 1125 1175.42857 NA NA Texas State
2020-05-04 1114 1167.14286 NA NA Texas State
2020-05-05 1258 1185.57143 NA NA Texas State
2020-05-06 1322 1191.85714 NA NA Texas State
2020-05-07 1089 1183.14286 NA NA Texas State
2020-05-08 1327 1189.42857 NA NA Texas State
2020-05-09 1623 1265.42857 NA NA Texas State
2020-05-10 1333 1295.14286 NA NA Texas State
2020-05-11 1256 1315.42857 NA NA Texas State
2020-05-12 1303 1321.85714 NA NA Texas State
2020-05-13 1410 1334.42857 NA NA Texas State
2020-05-14 1056 1329.71429 NA NA Texas State
2020-05-15 1550 1361.57143 NA NA Texas State
2020-05-16 1245 1307.57143 NA NA Texas State
2020-05-17 1392 1316.00000 NA NA Texas State
2020-05-18 1219 1310.71429 NA NA Texas State
2020-05-19 1353 1317.85714 NA NA Texas State
2020-05-20 1613 1346.85714 NA NA Texas State
2020-05-21 1510 1411.71429 NA NA Texas State
2020-05-22 1446 1396.85714 NA NA Texas State
2020-05-23 1472 1429.28571 NA NA Texas State
2020-05-24 1652 1466.42857 NA NA Texas State
2020-05-25 1285 1475.85714 NA NA Texas State
2020-05-26 1253 1461.57143 NA NA Texas State
2020-05-27 1295 1416.14286 NA NA Texas State
2020-05-28 1497 1414.28571 NA NA Texas State
2020-05-29 1525 1425.57143 NA NA Texas State
2020-05-30 1671 1454.00000 NA NA Texas State
2020-05-31 1540 1438.00000 NA NA Texas State
2020-06-01 1302 1440.42857 NA NA Texas State
2020-06-02 1656 1498.00000 NA NA Texas State
2020-06-03 1230 1488.71429 NA NA Texas State
2020-06-04 1401 1475.00000 NA NA Texas State
2020-06-05 1510 1472.85714 NA NA Texas State
2020-06-06 1584 1460.42857 NA NA Texas State
2020-06-07 1348 1433.00000 NA NA Texas State
2020-06-08 1439 1452.57143 NA NA Texas State
2020-06-09 1827 1477.00000 NA NA Texas State
2020-06-10 1335 1492.00000 NA NA Texas State
2020-06-11 1533 1510.85714 NA NA Texas State
2020-06-12 1763 1547.00000 NA NA Texas State
2020-06-13 1707 1564.57143 NA NA Texas State
2020-06-14 1653 1608.14286 NA NA Texas State
2020-06-15 1812 1661.42857 NA NA Texas State
2020-06-16 1805 1658.28571 NA NA Texas State
2020-06-17 2071 1763.42857 NA NA Texas State
2020-06-18 2372 1883.28571 NA NA Texas State
2020-06-19 2220 1948.57143 NA NA Texas State
2020-06-20 2408 2048.71429 NA NA Texas State
2020-06-21 2361 2149.85714 NA NA Texas State
2020-06-22 2358 2227.85714 NA NA Texas State
2020-06-23 2336 2303.71429 NA NA Texas State
2020-06-24 2762 2402.42857 NA NA Texas State
2020-06-25 2615 2437.14286 NA NA Texas State
2020-06-26 3150 2570.00000 NA NA Texas State
2020-06-27 2918 2642.85714 NA NA Texas State
2020-06-28 3184 2760.42857 NA NA Texas State
2020-06-29 2711 2810.85714 NA NA Texas State
2020-06-30 3540 2982.85714 NA NA Texas State
2020-07-01 3502 3088.57143 NA NA Texas State
2020-07-02 3586 3227.28571 NA NA Texas State
2020-07-03 4149 3370.00000 NA NA Texas State
2020-07-04 4232 3557.71429 NA NA Texas State
2020-07-05 3829 3649.85714 NA NA Texas State
2020-07-06 4198 3862.28571 NA NA Texas State
2020-07-07 4462 3994.00000 NA NA Texas State
2020-07-08 4333 4112.71429 NA NA Texas State
2020-07-09 4953 4308.00000 NA NA Texas State
2020-07-10 5308 4473.57143 NA NA Texas State
2020-07-11 4490 4510.42857 NA NA Texas State
2020-07-12 5408 4736.00000 NA NA Texas State
2020-07-13 6211 5023.57143 NA NA Texas State
2020-07-14 6117 5260.00000 NA NA Texas State
2020-07-15 4810 5328.14286 NA NA Texas State
2020-07-16 5520 5409.14286 NA NA Texas State
2020-07-17 4881 5348.14286 NA NA Texas State
2020-07-18 NA 5398.91800 5356.999 5440.837 Texas State
2020-07-19 NA 5449.69314 5377.357 5522.029 Texas State
2020-07-20 NA 5500.46828 5395.558 5605.378 Texas State
2020-07-21 NA 5551.24342 5411.035 5691.451 Texas State
2020-07-22 NA 5602.01856 5423.779 5780.258 Texas State
2020-07-23 NA 5652.79370 5433.889 5871.698 Texas State
2020-07-24 NA 5703.56884 5441.484 5965.653 Texas State
2020-07-25 NA 5754.34398 5446.681 6062.007 Texas State
2020-07-26 NA 5805.11912 5449.584 6160.654 Texas State
2020-07-27 NA 5855.89426 5450.292 6261.496 Texas State
forecast.plot(texas.arima.df.1, "Texas Daily Hospitalizations", "Daily Hosp")
## Warning: Removed 1 rows containing missing values (geom_path).

## Warning: Removed 136 rows containing missing values (geom_path).

## Warning: Removed 136 rows containing missing values (geom_path).

STANDARD STATISTICAL TESTS

CASE RATIOS

TSA

#TODO: update all varnames
cumcases.TSA <- TSA.Daily %>% dplyr::select(Date,Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.TSA$Date <- ymd(cumcases.TSA$Date)

latestdate <- max(cumcases.TSA$Date)
earliestdate <- latestdate - 14
middate <- latestdate - 7

week2 <- subset(cumcases.TSA, Date==latestdate, select=Cases_Cumulative) - 
  subset(cumcases.TSA, Date==middate, select=Cases_Cumulative)

week1 <- subset(cumcases.TSA, Date==middate, select=Cases_Cumulative) - 
  subset(cumcases.TSA, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
#current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA
tsa_current_ratio_dat <- data.frame(TSA=subset(cumcases.TSA, Date==latestdate, select=Level_Name)[,1],
                                    current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

# add TSA name
tsa_current_ratio_dat <- merge(tsa_current_ratio_dat, unique(TSA.daily[, c('TSA', 'TSA_Name')]), by = 'TSA')

# combine TSA name
tsa_current_ratio_dat$TSA <- paste0(tsa_current_ratio_dat$TSA, ' - ', tsa_current_ratio_dat$TSA_Name)
tsa_current_ratio_dat$TSA_Name <- NULL

PHR

cumcases.PHR <- PHR.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.PHR$Date <- ymd(cumcases.PHR$Date)

latestdate <- max(cumcases.PHR$Date)
earliestdate <- latestdate - 14
middate <- latestdate - 7

week2 <- subset(cumcases.PHR, Date==latestdate, select=Cases_Cumulative) - 
  subset(cumcases.PHR, Date==middate, select=Cases_Cumulative)

week1 <- subset(cumcases.PHR, Date==middate, select=Cases_Cumulative) - 
  subset(cumcases.PHR, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
#current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA
phr_current_ratio_dat <- data.frame(Level_Name=subset(cumcases.PHR, Date==latestdate, select=Level_Name)[,1],
                                    current_ratio=current_ratio, current_ratio_cat = cut(current_ratio, breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

# add PHR name
phr_current_ratio_dat <- merge(phr_current_ratio_dat, unique(PHR.Daily[, c('PHR', 'Level_Name')]), by = 'Level_Name')

# combine PHR name
phr_current_ratio_dat$Level_Name <- paste0(phr_current_ratio_dat$PHR, ' - ', phr_current_ratio_dat$Level_Name)
phr_current_ratio_dat$PHR <- NULL

colnames(phr_current_ratio_dat)[1] <- 'PHR'

County

# declare vals
cumcases.county <- county.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.county$Date <- ymd(cumcases.county$Date)

latestdate <- max(cumcases.county$Date, na.rm=T)
earliestdate <- latestdate - 14
middate <- latestdate-7

# calc cumulative case ratios from 2 weeks ago and last week
week2 <- subset(cumcases.county, Date==latestdate, select=Cases_Cumulative) - 
         subset(cumcases.county, Date==middate, select=Cases_Cumulative)

week1 <- subset(cumcases.county, Date==middate, select=Cases_Cumulative) - 
         subset(cumcases.county, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]

# add contingencies 
current_ratio[week1<10] <- NA
current_ratio[current_ratio<0] <- NA

# output
# TODO: investigate select vs subset(...)[, 'County']
county_current_ratio_dat <- 
  data.frame(County=subset(cumcases.county, Date==latestdate, select=Level_Name)[,1],
             current_ratio=current_ratio,
             current_ratio_cat = cut(current_ratio,
                                     breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

Metro

cumcases.metro <- metro.Daily %>% dplyr::select(Date, Level_Name, Cases_Cumulative) %>% filter(!is.na(Level_Name))
cumcases.metro$Date <- ymd(cumcases.metro$Date)

latestdate <- max(cumcases.metro$Date, na.rm=T)
earliestdate <- latestdate - 14
middate <- latestdate-7

week2 <- subset(cumcases.metro, Date==latestdate, select=Cases_Cumulative) - 
         subset(cumcases.metro, Date==middate, select=Cases_Cumulative)
week1 <- subset(cumcases.metro, Date==middate, select=Cases_Cumulative) - 
         subset(cumcases.metro, Date==earliestdate, select=Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
current_ratio[current_ratio<0] <- NA

metro_current_ratio_dat <-
  data.frame(Metro_Area=subset(cumcases.metro, Date==latestdate, select=Level_Name)[,1],
             current_ratio=current_ratio,
             current_ratio_cat = cut(current_ratio,
                                     breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

State

cumcases.state <- state.Daily %>% dplyr::select(Date, Cases_Cumulative)
cumcases.state$Date <- ymd(cumcases.state$Date)

latestdate <- max(cumcases.state$Date, na.rm = TRUE)
earliestdate <- latestdate - 14
middate <- latestdate-7

week2 <- subset(cumcases.state, Date == latestdate, select = Cases_Cumulative) - 
         subset(cumcases.state, Date == middate, select = Cases_Cumulative)

week1 <- subset(cumcases.state, Date == middate, select = Cases_Cumulative) - 
         subset(cumcases.state, Date == earliestdate, select = Cases_Cumulative)

current_ratio <- (week2/week1)[,1]
current_ratio[current_ratio<0] <- NA

state_current_ratio_dat <- 
  data.frame(current_ratio=current_ratio,
             current_ratio_cat = cut(current_ratio,
                                     breaks=c(0,0.5,.9,1.1,1.5,2,4,8,max(current_ratio))))

export

write.csv(state_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/state_case_ratio.csv', row.names = F)
write.csv(tsa_current_ratio_dat,'statistical-output/standard-stats/case-ratios/tsa_case_ratio.csv', row.names = F)
write.csv(county_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/county_case_ratio.csv', row.names = F)
write.csv(metro_current_ratio_dat,'statistical-output/standard-stats/case-ratios/metro_case_ratio.csv', row.names = F)
write.csv(phr_current_ratio_dat,'statistical-output/standard-stats/case-ratios/phr_case_ratio.csv', row.names = F)

write.csv(state_current_ratio_dat, 'tableau/case-ratios/state_case_ratio.csv', row.names = F)
write.csv(tsa_current_ratio_dat,'tableau/case-ratios/tsa_case_ratio.csv', row.names = F)
write.csv(county_current_ratio_dat, 'tableau/case-ratios/county_case_ratio.csv', row.names = F)
write.csv(metro_current_ratio_dat,'tableau/case-ratios/metro_case_ratio.csv', row.names = F)
write.csv(phr_current_ratio_dat,'tableau/case-ratios/phr_case_ratio.csv', row.names = F)

grouping

colnames(county_current_ratio_dat)[1] <- 'Level'
colnames(metro_current_ratio_dat)[1] <- 'Level'
colnames(tsa_current_ratio_dat)[1] <- 'Level'
colnames(phr_current_ratio_dat)[1] <- 'Level'
state_current_ratio_dat$Level <- 'Texas'


county_current_ratio_dat$Level_Type = 'County'
metro_current_ratio_dat$Level_Type = 'Metro'
tsa_current_ratio_dat$Level_Type = 'TSA'
phr_current_ratio_dat$Level_Type = 'PHR'
state_current_ratio_dat$Level_Type = 'State'

combined_current_ratio_dat <- rbind(county_current_ratio_dat, metro_current_ratio_dat,
                                   tsa_current_ratio_dat, state_current_ratio_dat, phr_current_ratio_dat)

write.csv(combined_current_ratio_dat, 'statistical-output/standard-stats/case-ratios/stacked_case_ratio.csv',
          row.names = FALSE)

# write.csv(combined_current_ratio_dat, 'tableau/stacked_case_ratio.csv',
#           row.names = FALSE)

% CHANGE

create.case.test <-   function(level, dat, region){
  # creates the % difference in cases and tests and smooth line with CIs 
  # level: either "TSA", "county", or "metro". Note that "county" won't work for many counties unless have enough cases.
  # dat: dataset (e.g. "county", "metro", "tsa")
  # region: the region within the dataset (county, metro region, or tsa)
  
  if(level == "TSA"){
    dat <- subset(dat, TSA==region)
  }
   if(level == "PHR"){
    dat <- subset(dat, PHR==region)
  }
  if(level == "county"){
    dat <- subset(dat, County == region)
  }
  if(level == "metro"){
    dat <- subset(dat, Metro_Area==region)
  }

  # restrict data to first test date (for % test increase)
  first.test.date <- ymd("2020-04-22")
  dat$Date <- ymd(dat$Date)
  dat$Date2 <- as.numeric(ymd(dat$Date))
  
  # get 14 day rolling avg for instances where initial cases or tests are 0
  dat$cases_avg14 <- rollmean(dat$Cases_Daily, k=14, fill = NA, align = 'right', na.rm = TRUE)
  dat$tests_avg14 <- rollmean(dat$Tests_Daily, k=14, fill = NA, align = 'right', na.rm = TRUE)
  
  # restrict new df to first.test.date forward  
  slopedata.tests <- data.frame(subset(dat, select=c(Date, Date2, Cases_Daily,
                                                     Tests_Daily, cases_avg14, tests_avg14))) %>%
    subset(ymd(Date) >= ymd(first.test.date))
  
  tests_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=Tests_Daily))
  cases_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=Cases_Daily))
  
  # numeric date for gam                         
  tmpdata <- data.frame(Date2=slopedata.tests$Date2)

  if (cases_start != 0 & tests_start != 0) {
    
    # Calc splines if cases and tests will produce non Inf values
    slopedata.tests$newtestsY <- 100*(as.vector(slopedata.tests$Tests_Daily / tests_start) - 1)
    slopedata.tests$newcasesY <- 100*(as.vector(slopedata.tests$Cases_Daily / cases_start) - 1)
    
    # fit and predict w/ gam for spline vals
    cases.gam <- gam(newcasesY~s(Date2,bs="cs", k=5), data=slopedata.tests)
    tests.gam <- gam(newtestsY~s(Date2,bs="cs",k=5), data=slopedata.tests)
    
    cases.fit <- predict(cases.gam, tmpdata, se.fit=TRUE)
    tests.fit <- predict(tests.gam, tmpdata, se.fit=TRUE)
    
    # Use 95% CI for splines
    tmp.df <- data.frame(Date = slopedata.tests$Date,
                         Data_Type = 'Raw',
                         cases_percentdiff = slopedata.tests$newcasesY,
                         tests_percentdiff = slopedata.tests$newtestsY,
                         cases_percentdiff_spline = cases.fit$fit, 
                         cases_percentdiff_spline_lower = cases.fit$fit-1.96*cases.fit$se, 
                         cases_percentdiff_spline_upper = cases.fit$fit+1.96*cases.fit$se,
                         tests_percentdiff_spline = tests.fit$fit,
                         tests_percentdiff_spline_lower = tests.fit$fit-1.96*tests.fit$se,
                         tests_percentdiff_spline_upper = tests.fit$fit+1.96*tests.fit$se)    
  } else {
    
    # only calc % increase for regions w/ sparse cases
    cases_avg_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=cases_avg14))
    tests_avg_start <- as.numeric(subset(slopedata.tests, Date==first.test.date, select=tests_avg14))

    slopedata.tests$newcasesY <- 100*(as.vector(slopedata.tests$cases_avg14 / cases_avg_start) - 1)
    slopedata.tests$newtestsY <- 100*(as.vector(slopedata.tests$tests_avg14 / tests_avg_start) - 1)
    
    tmp.df <- data.frame(Date = slopedata.tests$Date,
                         Data_Type = '14-Day Average',
                         cases_percentdiff = slopedata.tests$newcasesY,
                         tests_percentdiff = slopedata.tests$newtestsY,
                         cases_percentdiff_spline = rep(NA, times = nrow(slopedata.tests)),
                         cases_percentdiff_spline_lower = rep(NA, times = nrow(slopedata.tests)),
                         cases_percentdiff_spline_upper = rep(NA, times = nrow(slopedata.tests)),
                         tests_percentdiff_spline = rep(NA, times = nrow(slopedata.tests)),
                         tests_percentdiff_spline_lower = rep(NA, times = nrow(slopedata.tests)),
                         tests_percentdiff_spline_upper = rep(NA, times = nrow(slopedata.tests)))  
    }
  
  return(tmp.df)  
}

State

state <- readxl::read_xlsx('combined-datasets/state.xlsx', sheet=1)
state.case.test.df <- create.case.test(level="State", state, NA)
write.csv(state.case.test.df, 'statistical-output/standard-stats/pct-change/state_pct_change.csv',
          row.names = FALSE)

TSA

tsa <- read.csv('combined-datasets/tsa.csv')

# Create tsa-level data (can optimize this code in the future, but it runs pretty quickly already)
# TODO: convert to gapply -> rbindlist
tsalist <- unique(tsa$TSA)
tmp <- create.case.test(level="TSA", tsa,tsalist[1])
tsa.case.test.df <- data.frame(TSA=rep(tsalist[1], nrow(tmp)), 
                               TSA_Name=rep(unique(tsa$TSA_Name[1])),
                               tmp)

for(i in c(2:length(tsalist))){
  tmp <- create.case.test(level="TSA", tsa,tsalist[i])
  tsa.case.test.df <- rbind(tsa.case.test.df, data.frame(TSA=rep(tsalist[i], nrow(tmp)), 
                                                         TSA_Name=rep(unique(tsa$TSA_Name[i])),
                                                         tmp))
}


## combine tsa name
tsa.case.test.df$TSA = paste0(tsa.case.test.df$TSA, ' - ', tsa.case.test.df$TSA_Name)
tsa.case.test.df$TSA_Name = NULL

write.csv(tsa.case.test.df, 'statistical-output/standard-stats/pct-change/tsa_pct_change.csv',
          row.names = FALSE)

PHR

phr <- read.csv('combined-datasets/phr.csv')

# Create phr-level data (can optimize this code in the future, but it runs pretty quickly already)
# TODO: convert to gapply -> rbindlist
phrlist <- unique(phr$PHR)
tmp <- create.case.test(level="PHR", phr,phrlist[1])
phr.case.test.df <- data.frame(PHR=rep(phrlist[1], nrow(tmp)), 
                               PHR_Name=rep(unique(phr$PHR_Name[1])),
                               tmp)

for(i in c(2:length(phrlist))){
  tmp <- create.case.test(level="PHR", phr,phrlist[i])
  phr.case.test.df <- rbind(phr.case.test.df, data.frame(PHR=rep(phrlist[i], nrow(tmp)), 
                                                         PHR_Name=rep(unique(phr$PHR_Name[i])),
                                                         tmp))
}


## combine phr name
phr.case.test.df$PHR = paste0(phr.case.test.df$PHR, ' - ', phr.case.test.df$PHR_Name)
phr.case.test.df$PHR_Name = NULL

write.csv(phr.case.test.df, 'statistical-output/standard-stats/pct-change/phr_pct_change.csv',
          row.names = FALSE)

Metro

Metro <- read.csv('combined-datasets/metro.csv')

metrolist <- unique(Metro$Metro_Area)

tmp <- create.case.test(level="metro", Metro,metrolist[1])
metro.case.test.df <- data.frame(Metro_Area=rep(metrolist[1], nrow(tmp)), tmp)

tmp <- create.case.test(level="metro", Metro,metrolist[2])
metro.case.test.df <- rbind(metro.case.test.df, data.frame(Metro_Area=rep(metrolist[2], nrow(tmp)), tmp))

write.csv(metro.case.test.df, 'statistical-output/standard-stats/pct-change/metro_pct_change.csv',
          row.names = FALSE)

County

county <- read.csv('combined-datasets/county.csv')
countylist <- unique(county$County)

tmp <- create.case.test(level="county", county,countylist[1])
county.case.test.df <- data.frame(County=rep(countylist[1], nrow(tmp)), tmp)

for(i in 2:length(countylist)) {
  tmp <- create.case.test(level="county", county,countylist[i])
  county.case.test.df <- rbind(county.case.test.df, data.frame(County=rep(countylist[i], nrow(tmp)), tmp))
}

# assess missing vals
missing_vals <- sum(is.na(county.case.test.df$cases_percentdiff)) + sum(is.na(county.case.test.df$tests_percentdiff))
recorded_vals <- sum(!is.na(county.case.test.df$cases_percentdiff)) + sum(!is.na(county.case.test.df$tests_percentdiff))

# 11.88% w/ 2020-05-01
# 12.71% w/ 2020-04-22
# 12% w/ 2020-05-15
missing_vals / (missing_vals + recorded_vals)
## [1] 0.1035614
write.csv(county.case.test.df, 'statistical-output/standard-stats/pct-change/county_pct_change.csv',
          row.names = FALSE)

grouping

colnames(county.case.test.df)[1] <- 'Level'
colnames(metro.case.test.df)[1] <- 'Level'
colnames(tsa.case.test.df)[1] <- 'Level'
colnames(phr.case.test.df)[1] <- 'Level'
state.case.test.df$Level <- 'Texas'

county.case.test.df$Level_Type = 'County'
metro.case.test.df$Level_Type = 'Metro'
tsa.case.test.df$Level_Type = 'TSA'
phr.case.test.df$Level_Type = 'PHR'
state.case.test.df$Level_Type = 'State'

combined.case.test.df <- rbind(county.case.test.df, metro.case.test.df,
                               tsa.case.test.df, state.case.test.df, phr.case.test.df)
write.csv(combined.case.test.df, 'statistical-output/standard-stats/pct-change/stacked_pct_change.csv',
          row.names = FALSE)

write.csv(combined.case.test.df, 'tableau/stacked_pct_change.csv',
          row.names = FALSE)

STACK COMBINATIONS

combined_current_ratio_dat$Date = max(combined.case.test.df$Date)
colnames(combined.arima.df.1)[5:6] = c('TS_CI_Lower', 'TS_CI_Upper')
colnames(combined.rt.df)[4:5] = c('RT_CI_Lower', 'RT_CI_Upper')

stacked_all = Reduce(function(x, y) merge(x, y, by = c('Level_Type', 'Level', 'Date'), all=TRUE),
       list(combined.case.test.df, combined_current_ratio_dat, combined.arima.df.1, combined.rt.df))

write.csv(stacked_all, 'tableau/stacked_all.csv', row.names = FALSE)

rolling average

Omitted fow now, Dr. Yaseen is processing this in Tableau